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Generation  of  Vertically  Incident  Seismograms 

By  L.  NEIL  FRAZER,  DAVID  L.  BATES,  and  ALBERT  J.  RUDMAN 


Abstract 

An  algorithm  for  synthetic  seismograms 
following  the  method  of  Kennett  (1981)  has 
been  implemented  in  FORTRAN  77.  Compu¬ 
tations  proceed  in  the  frequency  domain  and 
then  are  Fourier  invert^.  Reflections  of 
normally  incident  waves  and  their  multiples 
are  plotted  on  a  Versatec  plotter.  The 
algorithm  allows  easy  adjustment  to  other 
plotting  protocols.  A  brief  outline  of  the 
theory  and  a  program  listing  with  examples 
are  included. 

Introduction 

This  report  concerns  the  generation  of 
synthetic  seismograms  for  vertically  incident 
waves  on  a  multilayered  horizontal  medium. 
Superposition  of  plane  waves  allows  the 
simulation  of  a  localized  source.  Compute- 

Theory 

In  wave  propagation  in  the  presence  of 
absorption,  solution  to  the  wave  equation 
leads  to  a  complex  velocity  V(w)  that 

A  = 


tions  proceed  in  the  frequency  domain  to 
permit  inclusion  of  attenuation  by  a  complex 
velocity.  The  plane-wave  components  are 
related  to  reflection  and  transmission  within 
the  layers  by  a  standard  filter  theory 
approach  following  the  method  of  Kennett 
(1981).  An  expansion  of  the  response  into  a 
power  series  permits  the  user  to  choose  a 
complete-response  seismogram  or  to  designate 
the  number  of  multiples  to  be  used.  The  final 
seismogram  is  constructed  by  Fourier  inver¬ 
sion. 

The  algorithm  presented  here.  Program 
VISP,  is  modified  from  Kennett  for  the  case 
of  horizontally  layered  media  bounded  on  top 
and  bottom  by  half  spaces.  This  code  does 
not  include  a  free  surfrice,  and  therefore  the 
ringing  effect  of  surface  multiples  is  not 
included. 


depends  on  the  Q  of  the  medium.  The 
absorption  model  is  an  exponential  decay 
with  distance 


where  a  is  related  to  frequency  w  and  Q 

[wl 

a  »  - 

2cQ 

(2) 

and  c  is  the  phase  velocity. 

A  plane  compressional  wave  normally  of  displacement  and  stress  determines  the 

incident  on  two  semi-infinite  layers  generates  amplitudes  of  the  waves  in  terms  of  the  re- 

a  reflected  and  transmitted  wave.  Continuity  flection  and  transmission  coefficients  R  and  T 


PjVj  -  ^  ^ 

PjVj  ♦  Pj^j  *  Pj+l'j+l 
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where  pj  is  the  density  of  layer  j  and  vj  is  the 
compressional  velocity. 

For  n  layers  sandwiched  between  two  half 


spaces,  a  matrix  representation  of  the  upgoing 
and  downgoing  wave  Reids  U  and  D  is  given 
by 


(::) 


=  Qj  Q2 


’■  C::) 


(4) 


Q  is  a  wave  propagator  related  to  a  matrix  of 
the  total  reflection  and  transmission  coeffi¬ 


cients  R  and  T  (Kennett,  1981,  equation 
3.81) 


Q  = 


Au  “ 

V'^DRu 

Td-i  - 

(5) 


Figure  1.  Graphic  representation  showing  the  first  few 
terms  of  equation  9  for  reflection  and  transmission 
matrices  of  superposed  media  and  showing  schemat¬ 
ically  the  interactions  undergone  by  the  waves  with 
the  regions  Zi  <  z  <  Zj  and  Z3  <  z  <  Zj.  From 
Kennett  (1981). 


where  the  subscripts  U  and  D  refer  to 
computations  for  upgoing  and  downgoing 
waves. 

To  illustrate  the  significance  of  the  total 
reflection  and  transmission  coefficient,  first 
consider  two  layers  sandwiched  between  half 
spaces  (fig.  1).  Then  the  overall  response  in 
terms  of  reflection  and  transmission  proper¬ 
ties  is 


p  13  —  12  A  1 2 p  23  fi  ^  ^2p  23"i“l  ip  12 

“d  ■  LI  ■  %  j  Id 

Ti  13  —  'p  23  Ft  p12p23t“1t>12 

Td  -  Tq  LI  -  %  I^D  J  Td 

Ru^3  =  r^23  ^  T23Ru^2  [j  . -r^2 35-^123 -1  ^^23 

Tu^^  =  *:pu^2  -  ■Rjj2SRyI23 -1  Yy2  3 


(6) 


n  o  o 
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COMMON  /  CHAR  /  IT  I  ME , I  DATE , lOBT , T 1 TLE 

W'H1TE(  6,  ’  (  IX,  "VISP  OCEAN  CRUST  MODELING  PROGRAM"  ,  8 X, 

»  A10,2X,A10///)' )  ITIME.IDATE 

WRITE(6,'("  EPS  =  ",F6.4,"  SIGMA  = " , F5 . 2 , 3X, "YS 1 DE= " , F5 . 2 , 
*3X, "TENTH  =  ", F5 . 2, /)’ )EPS, SIGMA, YSIDE, TENTH 
WRITE(6, ' ( "  LAYERU", IX, "VELOCITY" , 4X, "Q",5X, "DENS ITY" , 1 X, 

+ "THICKNESS " , 7X, "DEPTH" , / , 9X, "KM/SEC" , 2X, "KM/SEC" , 3X, "GRM/SEC" , 
-5X,"M",13X,"M")’ ) 

WRITE(6,'("  1  ",F8. 2,F9. l,F8.2) ' )  CL( 1 ) ,QL( 1 ) ,RHO( 1 ) 

TTOT  =  0.0 
DO  40  I=2,NL-1 

V\RITE(6,  ’  (  IX,  14  .  1X,F8.2  ,F9. 1  ,F8.2, 1X,F9.2  ,F10 . 2  , 

-  "  TO", 1X,F7.2) ' ) I ,CL( I ) ,QL( I ) ,RHO( I ) ,T( I ) ,TTOT,TTOT+T( I ) 

TTOT  =  TTOT+T(I) 

4  0  COST  I NUE 

I\RITE(6,  ’  (IX,  I4,1X,F8.2,F9.1,F8.2,F2  0.2, 

-  "  TO  INFINITY")')  NL,CL(NL) ,QL(NL) ,RHO(NL) ,TTOT 
RETURN 

END 


SL’BROUTINE  PLTARMI 

C 

(; - 


SL'B ROUTINE  PLTARMI 

COMMON  /  NUM  /  CL ( 2 0 0 ) , QL { 2 00 ) , RHO( 2 0 0 ) , T( 2 00 ) , EPS , S I GMA , 

*  CRUST , NL , NTRACE , LOBS , NW, TSEC , TSC , ASC , NMULT , 

*  TLAG,YS IDE, TENTH 
COMMON  /  CHAR  /  ITIME, IDATE, lOBT, TITLE 

COMMON  /FILTDT/  FI ( 1 0 ) , I DB( 1 0 ) , I F I LTYP ( 1 0 ) , NF I LT , MPHASE , 

-  FFl(lO) 

COMMON  /AR/ lAGC, WINDOW 
LOGICAL  MPHASE 
CIIARACTER*80  TITLE 
CTLARACTER’IO  IT  IME ,  I  DATE  ,  L INE 
CH.ARACTER*!  lOBT 
NFILT  =  0 

READ( 1 , ’ (A2 ) ' )  WWW 

READ( 1 , • )  NTRACE, NW, TSEC, TSC, ASC, NMULT, lAGC, 

•WINDOW 
WRITE (6, 101) 

10  I  FORMAT! / / , IX, "NTRACE", I  X, "#  FREQS", IX, "T( MSECS ) ",  IX, "IN/MSECS", 

•  IX,  "AMP-SCL",  IX,  "MULTPL",  IX,  "A(jC",  IX,  "WINDOW") 

WRITE!  6  ,  102  )NTRACE,NW, TSEC, TSC, ASC,NT<1ULT,  lAGC,  WINDOW 

102  FORMAT! 3 X, I  2 , 4 X, I  4 , 3 X, F6 . 0 , 4X, F4 . 3 , 6 X, F4 . 1 , 4 X, 1 2 , 3 X, 

* I3,2X,F4.0) 

IF  (NTRACE  .GT.  10. 0)  CALL  TERMIN!27,0) 

IF  ((NW.NE.16)  .AND.  (NW.NE.32)  .AND.  (NW.NE.64)  .AND.  (NW.NE.128) 

•  .AND.  (NW.NE.256)  .AND.  (NW  .NE.  512))  CALL  TERMIN! 23 ,NW) 

IF  (TSEC  .LT.  l.O)  CALL  TERMINI  25 ,NINT(TSEC) ) 

IF  ( (TSC.EQ.O.O)  .OR.  (TSC»TSEC. GT. YS I DE ) )  THEN 
TSC  =  YS IDE /TSEC 


onooo 
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IF  (NINTLAY+I  .GrT.  NL)  CALL  TERMIN(4  ,MLNR) 

TINTLAY  =  RHO(I) 

IF  (TINTLAY  .LT.  0)  THEN 
TINTLAY  =  -TINTLAY-CRUST 
IF  (TINTLAY  .LT.  0)  CALL  TERMIN( 5 ,MLNR) 

END  IF 

CRUST  =  CRUST+TINTLAY 

INCRCL  =  (CL(I+NINTLAY)-CL(I-1))/(NINTLAY+1) 

INCRQL  =  (QL(I+NINTLAY)-^(I-1))/(NINTLAY+1) 

IF  (RHO(  I+NINTLAY)  .EQ. -1 )  RHO(  I+NINTLA^'Is  (CL(  I +NINTLAY)  +  1 . 5  ) /3 
INCRRHO  =  (RHO( I+NINTLAY)-RHO( I-l ) )/(NINTLAY+1 ) 

DO  25  J=1,NINTLAY 

CL(J+I-1)  =  CL(I-l)  +  J* INCRCL 
QL(J+I-1)  =  ^(I-l)  +  J*INCRQL 
RHO(J+I-l)  =  RHO(I-l)  +  J*INCRRHO 
T(J+I-1)  =  TINTLAY/NINTLAY 
CONTINUE 
I  =  I+NINTLAY 
ELSE 

IF  (RHO(I)  .EQ.  -1)  RHO(I)=(CL(I)+1.5)/3 
IF  (T(I)  .LT.  0)  THEN 
T(I)  =  -T( I) -CRUST 

IF  (T(I)  .LT.  0)  CALL  TERMINI 5 ,MLNR) 

ENDIF 

IF  (I  .LT.  NL)  THEN 
CRUST  =  CRUST+Td) 

ELSE 

T(I)  =  10000 
ENDIF 
I  =  I  +  l 
ENDIF 

IFd.LT.NDQO  TO  20 

IF(RHO(I).EQ. -1)  RHOd)  =  (CL(I  )  +  l .  5  ) /3 
DO  30  1=2. NL 

IF  ((CL(I).LE.O)  .OR.  ((X(I).GE.IO))  CALL  TERMINI  1 1 , I ) 


IF  ((CL(I).LE.O)  .OR.  ((X(I).GE.IO))  CALL  TERMINI  1 1 , I ) 

IF  (<M.(I)  .LT.  0)  CALL  TERMIN(12,I) 

IFI I (RHOI I ) . LE. 0 ) .AND. (RHO( I ) .NE. -1 ) ) .OR. (RHO( I ) .GE. 10 ) ) 
CALL  TERMINUS,  I) 

CONTINUE 

END 


SUBROUTINE  PRTMLDA 


SUBROUTINE  PRTMLDA 
CHARACTER* 80  TITLE 
CHARACTER*10  ITIME, IDATE 
CHAKACrrER^  1  lOBT 

OONMON  /  NUM  /  CL(  200  )  ,^(  200  )  ,RHO(  200  )  ,T(  200  ) ,  EPS  ,S  lOVIA, 

►  GRUST , NL , NTRACE , LOBS , NW,  TSEC , TSC , ASC , NMULT , 

f  TLAG.YSIDE.TLNTH 
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END 


SUBROUTINE  MDLINCH 


SUBROUTINE  MDLINCH 
CHARACTER* 80  TITLE 
CHARACTER*10  ITIME.IDATE 
CHARACTER*!  lOBT 

COIVMON  /  NUM  /  CL(200).QL(200),RHO(200).T(200  ), EPS, SIGMA, 

♦  CRUST , NL , NTRACE , LOBS , NW, TSEC , TSC , ASC , NMULT , 

♦  TLAG,YSIDE,TLNTH 

COMVDN  /  CHAR  /  ITIME, IDATE, lOBT, TITLE 
REAL  INCRCL,  INCR^,  INCRRHO,TINTLAY 

TINTLAY=0.0 

ILNR=0 

READ  IN  TITLE  AND  WRITE  IT 

READ(1,100)TITLE 
100  FORMAT(A) 

WRITE(6,200)TITLE 
200  FORMAT!///, IX, A) 

NL=-1 

ILNR  =  ILNR+1 
DO  in  1=1,200 
NL  =  NL+1 

IF  (NL  .EQ.  0)  THEN 
READ(1,'(A2)')WWW 
READd,*)  EPS,SIGMA,YSIDE,TLNTH 
READ(1,’(A2)')WWW 
ELSE 

READd,*)  CL(NL),QL(NL),RHO(NL),T(NL) 

IF(CL(NL) .EQ.9999)00  TO  15 

. . . INTERPOLATED  LAYERS 

IF  (CL(NL)  .EQ.  0)  NL=NL+NINT(QL(NL) ) -1 
IF  (NL  .GT.  200)  CALL  TERMIN( 1 , I LNR) 

END  IF 

10  CONTINUE 

5  IF  ((CLd)  .EQ.  0)  .OR.  (T(l)  .LT.  0))  CALL  TERMIN(3,1) 

NL  =  NL-1 
CRUST  =  0 
MLNR= 1 
1  =  2 

...EXPAND  INTERPOLATED  LAYERS  AND  CONVERT  ABSOLUTE  DEPTHS  TO  THICKNESS. 

0  MLNR  =  MLNR+1 

IF  (CL( I )  .EQ.  0)  THEN 
NINTLAY  =  NINT(QL( I ) ) 


oooo  o o o o o o o o o o o n o o o n o o o o  ooooooooooo 
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THE  OCHIE  HAS  A  MAIN  PROGRAM  THAT  CALLS  6  SUBROUTINES.  THESE 
SUBROUTINES,  IN  TURN,  CALL  6  OTHER  SUBROITTINES.  OUTPITT  OCMSISTS 
OF  A  BRIEF  LINE  PRINTER  SIAMARY  OF  THE  INPOT  PARAMETERS  AND  A 
PLOT  OF  THE  SYNTHETIC  SEISMOGRAM  AND  LAYER  PARAMETERS  (VELOCITY, 
Q  AND  IKNSITY). 

ACTION  OF  PROGRAM: 


THE  PROGRAM  READS  AND  CHECKS  THE  MODEL  PARAMETERS.  IN  ALL  CASES  IT  THEN 
PLOTS  THE  PARAMETER  VALUES  AS  THEY  VARY  WITH  DEPTH.  THEY  ARE  PLOTTED  ON 
THE  SAME  GRAPH,  WITH  VALUES  OF  Q  REDUCED  BY  1000.  THE  Y  AXIS  FILLS  THE 
PAPER  WIDTH  AND  HAS  VALUES  0  -  MAX  PARAMETER  VALUE.  THE  LENGTH  OF  THE  X 
AXIS  IS  DEPENDENT  ON  BOTH  THE  TOTAL  THICKNESS  AND  THE  NU>BER  OF  LAYERS 
REQUESTED  IN  THE  MODEL. 

THE  MODEL  USES  THE  THEORY  OF  B . L . N . KENNETT  (ADVANCES  IN  APPLIED 
MECHANICS,  VOLUME  21,  PP. 79-167,  ACADEMIC  PRESS,  1981)  TO  COMPUTE 
RESPCWSE  OF  MODEL  FOR  SELECTED  FREQUENCIES  (SUBROUTINE  EXEMODL) 

2*PI/TSEC . NW*2*PI/TSEC.  THE  RESPONSE  VALUES  ARE  THEN 

FOURIER  TRANSFORMED  TO  THE  TIME  DOMAIN  AND  PLOTTED. 

IF  NO  TIME  SCALE  FACTOR  IS  SPECIFIED  (OR  IF  THE  ONE  GIVEN  IS 
INAPPROPRIATE),  A  VALUE  OF  1  MSEC  PER  INCH  OR  LESS,  IF  NECESSARY,  IS  USED. 
TRACES  ARE  SPACED  0.65  INCHES  APART  AND  ARE  SCALED  TO  A  MAXIMUM  MAGNITUDE 
OF  0.5  INCHES,  UNLESS  AN  ASC  VALUE  IS  DESIGNATED  (OTHER  THAN  1).  IN  THAT 
CASE,  PEAKS  ARE  TRUNCATED  AT  0.5  INCHES. 


MAIN  PROGRAM 


PROGRAM  TEST ( OUTPUT, TAPE6=OUTPUT) 

CHARACTER* 80  TITLE 

CHARACTER*10  DATE, TIME, ITIME, IDATE 

CHARACTER* 1  lOBT 

OOMVON  /  NUM  /  CL(200),QL(200),RHO(200),T(200), EPS, SIGMA, 

♦  CRUST , NL , NTRACE , LOBS , NW, TSEC , TSC , ASC , NMULT , 

♦  TLAG,YSIDE,TLNTH 
OOMVDN  /  CHAR  / ITIME, IDATE, lOBT, TITLE 
COMVCN  /AR/ I AGC,  WINDOW 
OPEN(UNIT=l,STATUS='OLD' ,FILE='MDL' ) 

ITIME=TIME( ) 

IDATE=DATE( ) 

CALL  MDLINCH 
CALL  PRTMLDA 
CALL  PLTARMI 
CALL  EXEMODL 

IF(  lAGC.EQ.  DCALL  AGCFREQ 

CALL  WRTPLT  r 


u  u  u  u  u 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


AMPLIFY  ALL  EVENTS  AND  CREATE  A  RINGING  EFFECT.  TOO  LARGE 
A  WINDOW  EFFECTIVELY  CANCELS  THE  AGCACTION. 

8)  CAPTIONS  (NOT  READ  BY  THE  COMPUTER) 

9)  PARAMETERS  SELECTING  DETECTOR  AND  LAYER  TO  PUT  IT 

LINE, LOBS 

WHERE 

LINE  -  IS  EITHER  "REFLECTION",  "VERTICAL"  OR  "PRESSURE" 
DEPENDING  ON  RECEIVER  TYPE  DESIRED 

LOBS  -  IS  LAYER  NUNBER  TO  PLACE  RECEIVER  (NUNBER  MUST 
BE  IN  COL  28) 

IF  THE  RECEIVER  CHOICE  IS  THE  REFLECTION  COEFFICIENT 
ONLY  LAYER  1  MAY  BE  DESIGNATED.  THE  CALCULATION  IS 
MADE  AT  THE  BASE  OF  THE  HALF  SPACE.  IF  THE  RECEIVER  IS 
DESIGNATED  AS  VERTICAL(DISPLACEMENT)  OR  PRESSURE,  THE 
RECEIVER  IS  AN  OBS  OR  BOREHOLE  INSTRUMENT  AND  MUST  BE 
PLACED  IN  LAYER  2  OR  GREATER  (NOT  LAYER  1).  THE  OBS 
IS  AT  THE  TOP  OF  THE  LAYER  SPECI F I ED,BUT  WITHIN  THE 
LAYER.  FOR  EXAMPLE,  AN  OBS  IN  LAYER  2  IS  AT  THE  TOP 
OF  LAYER  2,  ESSENTIALLY  ACROSS  THE  INTERFACE  FROM  THE 
STANDARD  REFLECTION  COEFFICIENT  CALCULATION  IN 
LAYER  1. 

10)  CAPTIONS  (NOT  READ  BY  COMPUTER) 

11)  IF  THE  PLOT  IS  IS  TO  BE  FILTERED,  THE  FILTER  PHASE  IS  GIVEN 

I  PHASE 

I  PHASE  -  IF  "0"  THE  FOLLOWING  FILTER(S)  IS  ZERO  PHASE, 

IF  "1"  THE  FILTER(S)  ARE  MINIMUM  PHASE 

12)  SUBSEQUENT  DATA  LINES  (UP  TO  10)  ARE  FILTER  DEFINITIONS 

IFILTYP, IDB,FFl 

WHERE: 

IFILTYP  -  1  =  HIGH  CUT  (LOW  PASS)  BUTTERWORTH 

2  =  LOW  CUT  (HIGH  PASS)  BUTTERWORTH 

3  =  NOTCH  FILTER 

IDB  -  DECIBEL  REDUCTION  WHERE  POLES  =  DB/6  +  1 
FFl  -  CUT-OFF  FREQUENCY  IN  HZ 

IF  A  BAND  PASS  FILTER  IS  DESIRED,  TWO  DATA  LINES  ARE 
READ  IN  CONSECUTIVELY:  ONE  GIVING  THE  HIGH  CUT  VALUES 
THE  OTHER  THE  LOW-CUT. 


IF  A  NOTCH  FILTER  IS  DESIRED,  INPUT  IS  SIMILAR  TO 
BAND  PASS  (SEE  ABOVE),  EXCEPT  IFILTYP  =  3  FOR 
BOTH  INPUT  LINES 


USING  THE  PROGRAM: 


oooooooooooooooooooooooonooonoooooooooooooonooooooooooo 
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TO  ZERO.  THE  FINAL  LAYER  IS  THE  BOTTOM  HALF-SPACE  AND  IS  IDENTIFIED 
BY  A  THICKNESS  T  SET  EQUAL  TO  ZERO.  A  FLAG  OF  4  CONSECUTIVE  9999 'S 
CLOSES  OFF  THE  LAYER  PARAMETER  INPUT. 

FOR  EACH  LAYER.  OTHER  THAN  THE  FIRST, IF  RHO  IS  GIVEN  AS  A  -1,  A 
DEFAULT  VALUE  OF  (CL+1.5)/3  IS  TAKEN. 

EXTRAPCXATION  -  IT  MAY  BE  DESIRED  TO  HAVE  LAYER  PARAMETERS  CHANGE 
SMOOTHLY  FOR  A  CERTAIN  REGION.  TO  AVOID  HAVING  TO  HAND  CALCULATE  AND 
INPUT  VALUES  FOR  SUCH  A  REGION  DIVIDED  INTO  A  DISCRETE  NUMBER  OF 
LAYERS,  IT  IS  POSSIBLE  TO  RE^ST  THAT  BETWEEN  THE  PRIOR  AND  NEXT 
FULLY  DESCRIBED  LAYERS  THERE  WILL  BE  N  LAYERS  OF  EQUAL  THICKNESS  AND 
WITH  PARAMETERS  COMPUTED  BY  LINEAR  EXTRAPOLATION  BETWEEN  THE  VALUES 
FOR  THE  PRIOR  AND  NEXT  LAYERS.  THIS  OPTIW  IS  REQUESTED  BY  ENTERING 
A  LINE  OF  THE  FORMAT  O.N.T  WHERE  THE  0  IS  AN  ESCAPE  VALUE  FOR  THE 
NORMALLY  EXPECTED  CL,  N  IS  THE  NOS  OF  LAYERS  TO  BE  INSERTED,  T 
IS  THE  TOTAL  THICKNESS  OF  THOSE  LAYERS,  AND  THE  "1"  IS  A  DUMVIY. 

FOR  EXAMPLE, THE  VALUE  OF  CL  FOR  THE  2ND  INTERPOLATED  LAYER  WILL  BE: 

CL(PRIOR)  ♦  2*(CL(NEXT)  -  CL(PRIOR) ) / (N+ 1 ) 

AS  A  USEFUL  GUIDE,  IF  YOU  WANT  EACH  LAYER  1/4  WAVELENGTH  THICK  (AT 
THE  HIGHEST  FREQUENCY  TO  BE  USED),  CHOOSE  N  =  T*NW*4 / (CL'TSEC) 

ABSOLUTE  DEPTH  -  NORMALLY  T  SHOULD  BE  A  POSITIVE  VALUE 
SPECIFYING  THE  THICKNESS  OF  THE  LAYER.  IF,  HOWEVER,  IT  IS  PREFERtD 
TO  SPECIFY  THE  DEPTH  BELOW  THE  HALF  SPACX  TO  WHICH  THE  LAYER  WILL 
EXTEND  (ITS  STARTING  POINT  WILL  OF  (XIURSE  BE  I^*1EDIATELY  BELOW  THE 
PREVIOUSLY  DESCRIBED  LAYER),  A  NEGATIVE  VALUE  CAN  BE  GIVEN,  WHOSE 
MAGNITUDE  WILL  INDICATE  THAT  ABSOLUTE  DEPTH.  THIS  OPTION  CAN  ALSO  BE 
BE  USED  WITH  EXTRAPOLATED  LAYERS. 

6)  CVVPTIONS  (NOT  READ  BY  COMPUTER) 


7)  INPUT  PARAMETERS  FOR  COMPUTATIONAL  PURPOSES  IN  FREE  FORMAT  ARE: 
NTRACE,NW,TSEC,TSC,ASC,NMULT,  lACKJ, WINDOW 

WHERE: 

NTRACE  -NUNBER  OF  TRACES  PLOTTED  AS  DUPLICATES  ON  THE  SYNTHETIC 
NW  -  NOS  OF  FREQUENCIES  (il6,  POWER  OF  2)  TO  BE  USED  IN  FOURIER 
TRANSFC»tM  TO  TIME  DOMAIN. 

TSEC  -  MSECS  OF  TIME  SERIES  (STARTS  AT  -1,  -5  C«  -10  MSECS 
DEPENDING  ON  THE  VALUE  OF  TSEC).  TSEC  SHOULD  BE  MORE 
THAN  TIME  FOR  1-3  REFLECTIONS  FROM  BOTTOM  INTERFACE,  DEPENDING 
ON  WHETHER  YOU  ARE  REQUESTING  PRIMARIES  ONLY  OR  FULL  RESPUISE. 
TSC  -  TIME  SCALING  FACTOR( IN/MSEC) . IF  TSEC*TSC4YS IDE .THEN  IF 
TSC  IS  POSITIVE, A  VALUE  OF  TSC  IS  COMPUTED  =  YSIDE/TSEC,  WHILE 
IF  TSC  IS  NEGATIVE, PLOT  IS  TRUNCATED  AT  YS IDE /TSC  MILLISECCM4DS 
ASC-  AMPL  SCALE  FACTOR  (DEFAULT  =  il.O).  IF  i  1,  SIGNAL  IS 
AMPLIFIED,  BUT  THE  PLOT  IS  CLIPPED  AT  0.5  INCHES  AMPLITUDE. 
NMULT  -NOS  OF  MULTIPLE  REFLECTIONS  TO  BE  COMPUTED. 

VALUE  4  0  IMPLIES  FULL  RESPOISE  REQUIRED  (BY  MATRIX  INVERSION). 
lAGC-  SET  =  TO  1  IF  AUTOMATICQAIN  CX)NTROL  IS  DESIRED 
WINDOW  -  LENGTH  OF  AGCWINDOW  IN  MSECS.  A  SMALL  WINDOW 
(RELATIVE  TO  EXPECTED  TIME  INTERVALS  OF  REFLECTIONS)  WILL 
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Appendix  3.  FORTRAN  77.  Program  VISP  and  Program  PLTVISP 


c 

C  PROGRAM  VISP 

C 

C  (VERTICAL  INCIDENT  SEISMOGRAM  PROGRAM) 

C 

C  BY  L.N.  FRAZER,  D.L.  BATES  &  A.J.  RUDMAN  JULY  1983 

C 

C 

C 

C  USES  I/O  STREAMS  AS  FOLLOWS: 

C  . . 

C  1:  INPUT 

C  6:  OUTPUT,  INCLUDING  ERROR  MESSAGES 

C  4:  STORED  ON  FILE  PLTIN  FOR  PLOTTING  BY  PROGRAM  PLTVISP 

C 

C 

C  FORMATS  AND  DESCRIPTIONS  OF  INPUT  FILE: 

C  . . 

C  A)  DATATYPES  -  ALL  DATA  IS  ACCEPTED  IN  FREE  FORMAT,  WITH  ANY  EMPTY  FIELDS 
C  BEING  READ  AS  0. 

C 

C  B)  DATALINES  -  USEABLE  DATALINES  ARE  INTERSPERSED  WITH  CAPTIONS -LINES 

C 

C  1)  TITLE  (ALPHA-NUMERIC) 

C 

C  2)  CAPTIONS  (NOT  READ  BY  COMPUTER) 

C 

C  3)  EPS,SIGMA,YSIDE,TLNTH 

C 

C  WHERE: 

C 

C  EPS  -  PARAMETER  GOVERNS  VARIATION  OF  VELOCITY  WITH  FREQUENCY  W 

C  SIGMA  -  PARAMETER  GOVERNS  VARIATION  OF  VELOCITY  WITH  FREQUENCY  W 

C  YSIDE  -  MAXIMUM  INCHES  PERMITTED  FOR  LENGTH  OF  SYNTHETIC  PLOT 

C  TLNTH  -  FIXED  LENGTH  OF  PARAMETER  PLOT(IN  INCHES) 

C 

C  4)  CAPTIONS  (NOT  READ  BY  COMPUTER) 

C 

C  5)  CL.QL.RHO.T 

C 

C  WHERE; 

C 

C  CL  -  IS  COMPRESS lONAL  VELOCITY  OF  LAYER  IN  KM/SEC 

C  QL  -  IS  COMPRESS lONAL  Q  FACTOR 

C  RHO  -  IS  DENSITY  IN  GM/CC 

C  T  -  IS  THICKNESS  IN  METERS 

C 

C  ALL  SUBSEQUENT  DATA  LINES  ARE  READ  AS  DESCRIPTICWS  OF  CONSECUTIVE 

C  LAYERS  OF  THE  MODEL  (UP  TO  200).  THE  FIRST  LAYER  IS  THE  UPPER 

C  HALF-SPACE  (CORRESPONDING  TO  AN  OCEAN  LAYER)  WITH  THICKNESS  SET 
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RVRB2EW 

SC 

SIGMA 

SIGNI 

SQREFT 

SQREFW 

T 

TD 

TDON 

TDONP 

TDOR 

TEMP 

TEST 

TINTLAY 

TITLE 

TLAG 

TENTH 

TSC 

TSEC 

TTOT 

TU 

TUNO 

TUNOP 

UNIT 

W 

WARNED 

WIGH 

WINDOW 

WW 

WWW 

Y 

YSIDE 


RVRB2*EWTGH«TU 
FFT  variable 

Parameter  used  in  computing  variation  of  velocities  with  W 

FFT  variable 

Square  of  amplitude 

Frequence  equivalent  of  SQREFT 

Array  of  thickness  (in  meters),  one  for  each  layer 
Downward  transmission  coefficients  for  current  interface 
Downward  transmission  coefficients  from  ocean  bottom  to 
lower  half  space 
New  value  of  TDON 

Downward  transmission  coefficients  from  ocean  bottom  to 
receiver 

Temporary  variable  for  truncation 
Temporary  variable 

Thickness  of  each  interpolated  layer 
String  used  for  plotting 
Time  plots  will  begin  at  -TLAG  msec 
Fixed  length  of  parameter  plot  (in  inches) 

Time-scaling  factor  in  inches/msec 
Number  of  msec  of  time  series 
Total  thickness  (in  meters) 

Upward  transmission  matr ix  for  current  interface 
Upwards  transmission  matrix  from  lower  half  space 
New  value  of  TUNO 
CMPLX  (1.0,  0.0) 

Current  frequency  in  KHz 
Logical  variable 

Phase  lag  associated  with  propagation  through  vertical 
distance 

ACX:  window  length  in  msec 
AGC  filter 

Ruse  to  read  past  a  card 

Common  factor  used  in  routine  REFL 

Max  length  of  plot  (in  inches) 
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I  STEP 

ITIME 

ITW2 

IW 

J 

L 

LINE 

LOBS 

LX 

M 

\1AXABSV 

MD 

MPHASE 

MLNR 

MU 

NINTLAY 

NFILT 

NL 

NMULT 

NT 

NTRACE 

NW 

PI 

POLES 

PR 

PKESSFA 

QL 

R 

HOOEFF 

HU 

KDON 

KDONP 

REFT 

REFW 

KHO 

HliOl 

HH02 

RTBl 

RTB2 

RTB3 

HU 

RUNO 

RUNOP 

RURO 

RVHBl 

RVRBIEW 

RVRB2 


FFT  variable 
Current  time 

Number  of  points  in  AGC  filter 
Frequency  loop  counter 
Miscellaneous  loop  counter 
FFT  variable 

String  to  hold  lines  read  from  input  file  (allows  free 
format).  (See  lOBT. ) 

Layer  in  which  OBS  located  (1  if  no  OBS) 

FFT  variable 
FFT  variable 
F i 1  ter  var iable 

(See  Kennett,  1981,  p.  108,  equation  3.38.) 

Flag:  If  set.  Minimum-phase  filter  required; 

otherwise.  Zero  phase 

Counter  of  model  file  input  lines,  used  in  error  messages 
(See  Kennett,  1981,  p.  106,  equation  3.38.) 

Number  of  interpolated  layers  requested 
Number  of  filters  requested 

Total  number  of  layers  in  model  including  bounding  half 
spaces 

Number  of  multiple  reflections  to  be  calculated 
Number  of  time  steps  (=  NW*2) 

Number  of  traces  to  be  plotted  on  the  seismogram 
Number  of  frequencies  to  be  calculated 
3.  14  IS926S3 
Filter  variable 

Value  of  pressure  recorded  at  OBS 
Pressure  related  variable  for  given  layer 
Array  of  P-wave  Q  values,  one  for  each  layer 
Filter  variable 

Reflection  coefficient  or  OBS  response  for  current 
frequency  W  and  slowness 

Downward  reflection  matrix  for  current  interface 
Downward  reflection  matrix  from  ocean  bottom  to  lower 
half  space 
New  value  of  RDON 

Reflection  coefficient  (or  OBS  response)  in  time  domain 

Reflection  coefficient  (or  OBS  response)  in  frequency 

doma i n 

Array  of  densities,  one  for  each  layer 

Density  of  layer  above  the  current  interface 

Density  of  layer  below  the  current  interface 

Temporary  variable 

Temporary  variable 

Variable  in  computation  of  RCX>EFF 

Upward  reflection  coefficients  for  current  interface 
Upward  reflection  coefficients  from  lower  half  space  to 
ocean  bottom 
New  value  of  KUNO 

Upward  reflection  coefficients  from  receiver  to  ocean 
bottom 

Effect  on  downward  waves  of  all  multiples  in  current  layer 
RVRBl*EWIGH*TDON 

Effect  on  upward  waves  of  all  multiples  in  current  layer 
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Appendix  2.  Glossary  of  Variables  Used  in  Program  VISP 

(Variables  in  the  plotting  Program  PLTV ISP  are  not  listed  here) 


A 

AA 

ALl 

AL2 

ASC 

B 

BB 

BEX 

CO 

Cl 

CL 

CRUST 

CTEMP 

CW 

CX 

CURRMIN 

DELFREQ 

DELTAT 

DET 

EIW 

EPS 

EWIGH 

EWRD 

EWRUNO 

FI 

FFl 

FACTOR 

FILT 

FlLTl 

FRSTTIM 

GL 

GLl 


1 

lAGC 
I  DATE 
IDB 

IFILTYP 

IL 

INCRCL 

INCRQL 

INCRRHO 

INDEXl 

1NDEX2 

lOBT 


1  PHASE 


Dumny  argument  in  multiple  computation 
FFT  variable 

P-wave  slowness  of  the  layer  above  the  current  interface 
P-wave  slowness  of  the  layer  below  the  current  interface 
An^litude  scaling  factor  for  traces 
Dumny  argument  in  multiple  computation 
AGC  filter  variable 

Numbers  smaller  than  exp  (-BEX)  taken  to  be  zero 
CM>LX  (0.0. 0.0) 

CMPLX  (1.0. 0.0) 

Array  of  P-wave  velocities,  one  per  layer 

Total  thickness  of  layers,  excluding  first  and  last 

FFT  variable 

FFT  variable 

FFT  variable 

AGC  filter  variable 

Frequency  increment  in  filtering 

Time  increment 

Temporary  variable  in  multiple  computation 
CMPLX(EPS,-W)**S1GMA 

Parameter  governing  variation  of  velocities  with 
frequency  W 
EXP(WIGH) 

EWIGH'RD 
EWIGH* RUNO 

Cut  off  frequency  in  Khz 

Cut  off  frequency  in  Hz 

Filter  variable 

Digitized  filter 

Second  copy  of  filter 

Logical  variable 

Z-comp  of  P-slowness  vector 

Z-comp  of  P-slowness  vector  in  layer  above  current 
interface 

Z-comp  of  P-slowness  vector  in  layer  below  current 
interface 

Miscellaneous  loop  counter 
Flag  for  AGC  f i I  ter 
Current  date 

DB  reduction  for  notch  filter,  or  =6*POLES-l  for 
Butterworth 

Type  of  each  filter  requested 
Layer  loop  counter 

Incremental  value  of  CL  through  interpolated  layers 
Incremental  value  of  QL  through  interpolated  layers 
Incremental  value  of  RHO  through  interpolated  layers 
Filter  variable 
Filter  variable  ' 

Character  for  type  of  OBS  response:  reflection/ 
vertical/pressure 

Flag:  0  for  zero  phase  and  I  for  minimum  phase 
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i^pendix  1.  Generalized  Description  Diagram  of  Programs  VISP  and  PLTVISP 


Subroutine 


1.  MAIN 


2.  MDLINCH 


3.  PRTMLDA 


4.  PLTARMI 


5.  EXEMODL 


6.  AGCFREQ 


7.  WRTPLT 


Subroutine 


1.  MAIN 


2.  READINP 


3,  INITPLT 


4.  PARMPLT 


5.  RSPMDL 


Program  VISP 


Comments 


Calls  six  subroutines  listed  below. 


Reads  model  parameters  and  interpolates  layer  parameters.  Calls 
TERMIN  for  error  stop. 


Writes  model  parameters. 

Reads  and  writes  computational  parameters,  user  options,  and  niter 
specitications.  Calls  TERMIN  for  error  stop. 

Main  computation  subroutine  for  a  given  frequency:  Computes 
reflection  coefficient  for  all  layers  (for  user-designated  number  of 
multiples).  Options  for  pressure  or  vertical  displacement.  Repeats  for 
all  frequencies  and  then  filters.  Finally,  transforms  to  time  domain. 

The  above  computations  involve  calls  to  flve  subroutines:  PARMGEN 
(computes  vertical  slowness),  REFL  (computes  reflection  coefflcient 
for  a  single  interface  at  a  single  frequency),  REVERB  (selects 
multiple  computations),  FILTER  (filters  amplitudes),  and  FFT  (fast 
Fourier  transform). 

Simulates  automatic  gain  control  in  frequency  domain  for  specified 
window.  Calls  FFT. 


Writes  header  and  seismogram  amplitudes  on  file  for  plotting  by 
lYogtam  PLTVISP.  (See  table  below.) 


Program  PLTVISP 


Comments 


Calls  four  subroutines  listed  below. 


Reads  data  generated  by  Program  VISP. 

Plots  major  headings. 

Initializes  plot  parameters  and  plots  axes  and  model  parameters. 

Writes  filter,  AGC,  and  multiple  specifications.  Plots  seismogram  traces 


and  axes. 


'•  ’■  •  • '  *'  •  •**  »*r  V  •”**  -  .•  •  •  '*'•  **’*^»*^  -  *•  •'*  •  •  -  •  ►*“  •**  •  f”'  •**,>* 


LITERATURE  CITED 

typical  field  records.  As  a  demonstration^  a 
model  that  simulates  a  continuous  velocity 
log  with  123  layers  was  used  to  create  Test  1 
(fig.  2).  A  Q)C  Cyber  170/855  computer  at 
the  Indiana  University  Wrubel  Computing 
Center  generated  the  output  data  in  13.6 
seconds  with  a  71300  octal  field  length.  Table 
1  lists  the  input  data  used  to  generate  figure 
2. 

Appendix  4  lists  the  input  data  and  output 
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plots  for  12  other  test  cases  used  to 
demonstrate  the  user  options. 

Literature  Cited 

Kennett,  B.  L.  N. 

1981  -  Elastic  wave  propogstion  in  stratified 
media,  in  Advances  in  applied  mechanics: 
New  York,  Academic  Press,  Inc.,  v.  21,  p. 
79-167. 


Table  1.  Input  data  used  to  generate  a  typical  synthetic  seismogram  (fig.  2) 


TEST  1 

-  DSYN. 

SYNTHETIC 

SEISMOGRAM 

EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

12.0 

CL 

RHO 

THICKNESS 

1.5 

2000. 

1.1 

0. 

1.9 

22. 

2.2 

10. 

2.1 

22. 

2.2 

10 

2.0 

22. 

2.2 

10 

2.2 

22. 

2.2 

10 

2.1 

22. 

2.2 

10 

2.8 

25. 

2.5 

10 

2.9 

e 

25. 

e 

2.5 

e 

10 

• 

e 

e 

e 

e 

e 

e 

• 

e 

e 

e 

e 

• 

2.7 

23 

2.4 

10 

2.6 

23 

2.4 

10 

2.7 

23 

2.4 

10 

2.8 

23 

2.4 

10 

2.6 

23 

2.4 

10 

1.5 

20 

2.2 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  2000.  .006  1.0  -1  1  400 

DETECTOR  TYPE  DETECTOR  LOCATICW 
REFLECTICHf  LAYER  1 

FILTER  PHASE 
1  (MINIMUM  PHASE) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 

1  72  70 

2  72  30 
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VI5P:  PLOTS  OF  PARAMETERS 

TE5T  l  -  osw.  5tN^rtCTlC  SClSWAfl. 


TEST  I  -  OSrN.  SrNTMCTIC  SEISWGRAfl. 

RCFLCCTIOM  COEFnCICMT 

AT  085  IN  layer  I  at  .0  n  QCPTH 

FOEOrNHI-  Siz  MNSeCSI-ZOOO 

NIN  •  .so  HZ  HAX  >ZS6.30  HZ 

FXL  BE5P0H5E  lALL  fU-HPLESl 

AGC  HINOOM  •  AOO.rtSeCS 

--FILTERS  •  HiNirw^  phase - 

HIGH  CUT  AT  ro.  HZirz  oe  slope 

LON  CUT  AT  30.  hZjFZ  30  SLOPE 


B 


FIguR  2.  (Test  1)  Output  plots  generated  by  Program  VISP.  A,  input  data:  velocity  (CL), 
Q  (QL),  and  density  (RHO);  B,  synthetic  seisntognm  genented  hom  input  data. 
Signiflcant  parameten  an  listed  above  the  seismogram. 


SUMMARY 

vi^ere,  for  example,  the  coefficients 
for  the  interface  at  Z2,  and  the  coefficients 
are  for  the  interlace  at  Z3,  but  are 


3 

adjusted  for  travel  through  the  layer  by  E  > 
eiwd/c,  where  d  is  layer  thickness,  ^nce 


=  eRu^^e  _  er^^^e  =  et^^s  -7^23  =  ETy^^ 

(7) 


Kennett  (1981)  pointed  out  that  the  power  series 
bracketed  quantity  may  be  expanded  as  a 

[I  -  A]‘^  =  I  +  a  +  a2  +  .  .  . 


Therefore, 


(8) 

‘(9) 


which  is  represented  schematically  in  figure  1. 
RD^^  is  the  reflection  from  the  upper  region. 
Reading  from  right  to  left,  Tu“  RD*^  Td‘* 
represents  a  transmission  back  up  through  the 
upper  zone.  Tu‘*  Rd*^  Ru‘*  RD^^  Td‘* 
represents  a  reverberation  within  the  layer. 

The  total  response  (equation  6)  includes  all 
internal  reverberations,  and  the  bracketed 
quantity  is  termed  the  reverberation  operator. 
Note  that  the  reflection-transmission  coeffi¬ 
cient  may  be  calculated  iteratively  by  adding 
a  layer  at  each  stage,  beginning  at  the  top  or 
the  bottom  of  the  stack  of  layers.  Computa¬ 
tion  of  the  response  of  the  model  proceeds 
for  selected  frequencies  and  is  finally  Fourier 
transformed  to  the  time  domain  and  plotted. 

Summary 

The  algorithm  for  Program  VISP  and  Program 
PLTVISP  (appendix  1)  presumes  an  incident 
plane  wave  originating  at  the  base  of,  but 
within,  the  upper  half  space.  The  two 
programs  are  written  in  FORTRAN  77,  and 
the  codes  are  listed  in  appendix  3.  Only  one 
input  file  is  required  to  run  Program  VISP. 
One  output  file  is  generated  for  plotting. 
Program  PLTVISP  plots  the  output  on  a 
Versatec  plotter  by  using  standard  calcomp 
calls.  A  glossary  of  constants  and  variables  is 
in  appendix  2.  The  comment  cards  at  the 


beginning  of  Program  VISP  detail,  card  by 
card,  the  input  data  required  to  generate  an 
output. 

User’s  options  include  an  output  plot 
(seismogram)  of  the  reflection  coefficient 
(pressure  of  the  reflected  wave  divided  by  the 
pressure  of  the  incident  wave),  pressure,  or 
vertical  displacement.  The  receiver  may  be  at 
the  top  of  any  of  the  layers  except  when  a 
reflection  coefficient  is  requested;  then  the 
receiver  is  just  above  the  first  interface. 
Seismograms  may  be  filtered  by  using  a 
Butterworth-type  high-pass,  low-pass,  or 
band-pass  filter.  The  filter  slopes  and  the 
phase  characteristics  (minimum  or  zero)  are 
specified  in  the  input.  The  user  may 
uniformly  amplify  ^e  signal  (ASC)  or 
optionally  select  a  simulated  automatic  gain 
control  (AGC).  The  final  seismogram  may 
display  primary  reflections  only,  primaries 
and  a  selected  number  of  multiples,  or  a 
complete  response. 

Besides  computational  parameters,  the  user 
must  read  in  the  following  model  parameters 
for  each  layer  and  the  top  and  the  bottom 
half  space:  velocity,  Q,  density,  and  thickness. 
If  a  negative  one  is  read  in  for  the  density,  a 
default  value  of  (velocity  1.5)/3  is 
calculated.  Program  VISP  successfully  gener¬ 
ates  synthetic  seismograms  that  simulate 


APPENIHX3 

ELSE  IF  (TSC.LT.0.0)  THEN 

C 

C  ...IN  THIS  CASE  DONT  PLOT  ALL. 

C 

TSC  =  -TSC 

END  IF 

IF  (ASC  .EQ.  0.0)  ASC=1.0 

IF  ((ASC  .LT.  1.0)  .C«.  (ASC  .GT.  100))  CALL  TERMIN(29,0) 

C 

C  LOOK  FOR  OBS  OPTIONS 

C  . 

c 

READd  ,  ’  (A2) '  )WWW 

READ(1,'(A11,16X, Il)')LINE,LOBS 
WRITE(6, 103 )LINE. LOBS 

103  FORMAT (/, IX, "DETECTOR  TYPE  IS  ",A,”  IN  LAYER  #",I3) 
IOBT=LINE(l !l) 

IF  (( lOBT.EQ. 'V )  .OR.  ( lOBT.EQ. ’P* ) )  THEN 

IF  ((LOBS.LT.2)  .OR.  ( LC»S .GT.NL) )  CALL  TERMIN(31,0) 
ELSE 
END  IF 

C 

C  NOW  READ  FILTER  PARAMETERS 

C  . . 

c 

READd,  ’  (A2) '  )WWW 
READd,*)  IPHASE 
IF( IPHASE. EQ.0)THEN 

WRITE! 6, ' (IX, "FILTER-PHASE") • ) 

WRITE!  6,  '  dX,  II,  "(ZERO  PHASE) ")')  IPHASE 
MPHASE=. FALSE. 

ELSE  I F ( 1 PHASE . EQ. 1 )THEN 

WRITE(6, ' (IX, "FILTER-PHASE") • ) 

WRITE(S, '(IX, II, "(MINIMUM  PHASE) ")') IPHASK 
MPHASEs.TRUE. 

ELSEIF( IPHASE.NE.O.AND. IPHASE.NE. 1 )THEN 
WRITE(6, ' (IX, "NO  FILTERS")') 

GO  TO  15 
ELSE 
END  IF 

WRITE(6, '(//, IX, "FILTER  TYPE",2X,"DB  SLOPE", 2X, 

+  "CUT-OFF  FREQUENCY")') 

READd, '(A2)')WWW 
DO  10  1=1,15 

READ(1,*,END=15)  IFILTYP( I ) , IDB( I ) ,FF1 ( I ) 

IF( IFILTYP( I ) .EQ. 1 )aO  TO  11 
IF(IFILTYP(I).EQ.3)GO  TO  41 
00  TO  21 

11  WRITE(6,104)IFILTYP(I),IDB(I),FF1(I) 

104  FORMAT (IX, II, "(LOW  PASS ) " , 4X, 1 3 , lOX, F4 . 0 ) 

GO  TO  31 

21  WRITE(6,105)IFILTYP(1),IDB(I),FF1(I) 

105  FORMAT(lX,Il,"(HIGH  PASS ) " , 3X, 1 3 , 1  OX, F4 . 0 ) 

GO  TO  31 

41  WRITE(6,106)IFILTYP(I),IDB(I),FF1(I) 


non  oooo  oonoo 
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106  FORMATdX,  II,  "(NOTCH  FILTER)  " ,  3X,  1 3  ,  lOX,  F4 . 0 , 3X,  F4 . 0  ) 

31  NFILT=NFILT+1 

Fl(NFILT)  =  FFl(NFILT)*.001 
10  CONTINUE 

15  IF  (NFILT  .GT.  10)  CALL  TERMIN(30,NFILT) 

CLOSE{UNlT=l 

RETURN 

END 


SUBROUTINE  EXEMODL 


SUBROUTINE  EXEMODL 

COMMON  /  NUM  /  CL( 200 ) ,QL( 200 ) , RHO( 200 ) ,T( 200 ) , EPS ,S IGMA, 

♦  CRUST , NL , NTRACE , LOBS , NW, TS  EC , TSC , ASC , NMULT , 

+  TLAG,YSIDE,TLNTH 

COMVION  /  CHAR  /  ITIME,  IDATE,  lOBT, TITLE 

COMMON  /FILTDT/  Fl ( 10 ) , I DB( 10 ) , IF ILTYP ( 10 ) , NFILT.MPHASE , 

+  FFKlO) 

CHARACTER*!  lOBT 
CHARACTER‘10  ITIME, IDATE 
CHARACTER*80  TITLE 
COMPLEX  GL1,GL2 

COMPLEX  EIW,WIGH,ROOEFF,CO,Cl ,AL1 ,AL2 
COMPLEX  EWIGH.RVRBl ,RVRB2 
COMPLEX  EWRUN0,EVvRD,RVRB1EW,RVRB2EW 
COMPLEX  RTBl,RTB2,RTB3,TrONP,TUNOP 
COMPLEX  RD0NP,RUN0P,RD0N,RUN0 
COMPLEX  TD0N,TUN0,RD,RU,TD,TU 
COMPLEX  MU,MD,RUR0,TD0R 
COMPLEX  PRESSFA, PR, UNIT, REFW(  1024  ) 

COMMW  /TS/  REFT(  1024) 

DATA  PI/3.141592653/ 

DATA  BEX/ -20.0/ 

DATA  CO /(O. 0,0.0)/,  Cl /(I. 0,0.0)/ 

INITIALIZE  BEGINNING  OF  PLOT 


NT  =  2*NW 

IF  (TSEC  .GT.  150.0)  THEN 
TLAG  =  10.0 

ELSE  IF  (TSEC  .GT.  90.0)  THEN 
TLAG  =5.0 
ELSE 

TLAG  =  1.0 
ENDIF 

UNIT  =  Cl 

LOOP  OVER  ALL  FREQUENCIES 


oono  ono  ooo  noon 
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C 

DO  10  IW=1,  NW-l 
W  =  1W*PI*2/TSEC 
EIW  =  CMPLX(EPS, -Wl^'SICaVIA 

ALl  =  CMPLX(1.0/CL(1),0.0)*(CMPLX(0.5/QL(1),0.0)/EIW+C1) 
RHOl  =  RHO(l) 

CALL  PARMGEN(AL1 ,GL1 ) 

LOOP  OVER  ALL  LAYERS 


DO  15  IL=1,  NL-1 

AL2  =  CMPLX(1.0/CL(IL+1),0.)  • 

(CMPLX(0.5/QL(1L+1),0.0)/EIW  +  Cl) 

RH02  =  RHO(IL+l) 

CALL  PARMGEN(AL2,GL2) 

CALL  REFL ( RHOl , RH02 , GLl , GL2 , RU , RD, TU , TD ) 

IF  ( I L  . EQ.  1 )  THEN 
RD0N=RD 
RUN0=RU 
TD0N=TD 
TUN0=TU 
MD  =  CO 
MU  =  CO 
ELSE 

WIGH  =  CMPLX(0.0,W*T( IL) )*GL1 
IF  (REAL(WIGH)  .GT.  BEX)  THEN 
EWIGH  =  CEXP(WIGH) 

ELSE 

EWIGH  =  CO 
END  IF 

EWRD=EWIGH«RD 
EWRUN0=EWiai*RUN0 
RTB1=EWRUN0*EWRD 
IF  (NMULT  .LT.  0)  THEN 

...COMPOTE  ALL  MULTIPLES  RVRB2=l+EWIGH*RD*RVRBl*EWIGH*RUNO. 

CALL  REVERB (RTBl.RVRBl) 

RTBl =RVRB1 • EWRUNO 
RTB2=EWRD*RTB1 
RVRB2=UNIT+RTB2 
ELSE 

...VALUE  FOR  NO  MULTIPLES. 

RVRB1=UNIT 

RVRB2=UNIT 

IF  (NMULT  .GT.  0)  THEN 
...COMPOTE  FINITE  NUMBER  OF  MULTIPLES. 


DO  20  1=1,  NMULT 
RTB2=RVRB1*RTB1 
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20 

RVRB1=UNIT+BTB2 
20  CONTINUE 

RTBl=EWRD*EWRUNO 
DO  25  1=1,  NMULT 
RTB2=RVRB2*RTB1 
RVRB2=UNIT+RTB2 
25  CONTINUE 

ENDIF 
ENDIF 

RTB1=EWIGH*TD0N 
RVRB1EW=RVRB1*RTBI 

C 

C  . . .TD0N=TD*RVRBl*EWIGH*TD0N. 

C 

TD0NP=TD*RVRB1EW 
RTB1=EWIGH*TU 
RVRB2  EW=  RVRB2  *  RTB I 

C 

C  . . .TUN0=TUN0*RVRB2*£WICM»TU 

C 

TUN0P=TUN0 •RVRB2EW 
RTB 1 = EWRD* RV RB 1 EW 

C 

C  . . .BD0N=RD0N+TUN0*EWIGH*RD*RVRB1*EWIGH*TD0N. 

C 

RTB2=TUN0*RTB1 
RD0NP=RD0N+RTB2 
RTB 1 = E WRUNO • RVRB 2EW 
C 

C  ... RUN0=RU+TD*EWIGH*RUN0*RVRB2*EWIGH*TU. 

C 

RTB2=TD*RTB1 

RUN0P=RU+RTB2 

TD0N=TD0NP 

TUN0=TUN0P 

RDON=RDONP 

RUN0=RUN0P 

ENDIF 

TEST  =  CABS(TDON) 

C 

C  . . .EXIT  LAYER  LOOP 

C 

IF((TEST.EQ.O.O).OR.(ALOGIO(TEST).LT.'8))GO  TO  30 
IF  (IL  .EQ.  LOBS-1)  THEN 
RL'R0  =  RUN0 
TD0R=TD0N 
TD0N=UNIT 
TUN0=UNIT 
RUNO  =  CO 
RUON  =  CO 

PRESSFA  =  RHO(LOBS)*(I.0-(4. 0/3.0)) 

C 

C  ...THE  OBS.  SITS  IN  THE  UPPERMOST  PART. 

C 

\tu  =  -GL2 


v--. 

•  -  • 


A 
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MD  =  GL2 
END  IF 

ALl  =  AL2 
RHOl  =  RH02 
GLl  =  GL2 

. . . END  OF  THE  LAYER  LOOP . 

15  CONTINUE 

30  CONTINUE 

IF  (LOBS  .GT.  1)  THEN 

...IF  AN  OBS  MOTION  IS  ASKED  FOR  THE  INCIDENT  PRESSURE  WAVE 
...IS  A  DERIVATIVE  OF  A  DELTA  FUNCTION.  IF  AN  OBS  PRESSURE 
...IS  ASKED  FOR,  THE  INCIDENT  PRESSURE  WAVE  IS  ASSUMED  TO  BE 
...A  DELTA  FUNCTION. 

RTB1=RUR0*RD0N 
IF(NMULT.LT.O)  THEN 

CALL  REVERB(RTB1 ,RTB2) 

ELSE 

RTB2=UNIT 

IF(NMULT.GT.O)  THEN 
DO  40  I=1,NMULT 

RTB2=UNIT*RTB1*RTB2 
)  CONTINUE 
ENDIF 
ENDIF 

RTB1=RTB2*TD0R 
RTB2=RD0N*RTB1 
PR=PRESSFA*(RTB2+RTB1 ) 

RTB2=MU*RD0N 

RTB2=MD+RTB2 

RTB3=RTB2*RTB1 

ADJUST  FOR  RECEIVER  TYPE 


IFdOBT.EQ.  'V  )THEN 
RCOEFF  =  RTB3/RHO(l) 

ELSE 

. . .PRESSURE. 

RCOEFF  =  PR/RHOd) 

ENDIF 

ELSE 

RCOEFF  =  RDON 
ENDIF 

RCOEFF  =  RCOEFF  •  CMPLX( 50 . • ( 1 +005 (P I* I W/NW) ) , 0 . )  • 
CEXP ( CMPLX( 0 . 0 , W*TLAG) ) 

REFWdW*!)  =  RCOEFF 

...  NEGATIVE  FRE<3UENCIES  ARE  COMPLEX  CONJ. 


nnonoonooo  oonoo  ooocd  ooc  non  onn 
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REFW(NT-IW*1 )  =  CONJG(RCOEPF) 

. . . END  FREQUENCY  LOOP 
0  CONTINUE 

. . .DC  AND  NYQUIST  FREQUENCY  SET  =  TO  0 

REFW(l)  =  CO 
REFW(NW*1)  =  REFW(l) 

CALL  FILTER(REFW,NW,TSEC) 

. . .TRANSFORM  TO  TIME  DOMAIN 

CALL  FFTINT.REFW, -I . ) 

DO  35  1=1 ,  NT 

REFT(I)  =  REAL(REFW( I ) ) 

5  CONTINUE 
RETURN 
END 


SUBROUTINE  PAR.VIGEN 


SUBROUTINE  PARMGEN ( AL , GL ) 
COMPLEX  AL,GL 
GL  =  CSQRTIAL'AL) 

IF  (AIMAG(GL).LT.O.O)  GL= -GL 

RETURN 

END 


SUBROUTINE  REEL 


SUBROUTINE  REFL( RHOl , RH02 ,GLl ,GL2 , RU, RD,TU,TD) 

...COMPUTES  REFLECTION  COEFFICIENT  RU,RD,TU  *  TD  FOR  GIVEN  LAYER 
...INTERFACE  PARAMETERS.  THE  REFLECT lON/TRANSMlSS ION  COEFFICIENTS  ARE 
...RATIOS  OF  DISPLACEMENT  POTENTIALS.  IN  THE  KTH  LAYER,  THE  DISPLACEMENT 
...VECTOR  U(K)  IS  GIVEN  BY: 

U(K)  =  GRAD(PHIU(K)«^PHID(K) )/(  I*W)  ^ 

CURL  CURL(ZHAT(PSIU(K)^PSID(K)))  /(W*W*P) 

...  IN  WHICH  PHIU  AND  PHID  ARE  THE  POTENTIALS  ASSOCIATED  WITH  THE  UP  AND 
...DOWN  GOING  P-WAVES  RESPECTIVELY.  ZHAT  IS  A  UNIT  VECTOR 
...POINTING  IN  THE  POSITIVE  Z  DIRECTION,  I.E.  DOWNWARDS. 

...THEN,  FOR  EXAMPLE,  THE  DOWNWARD  TRANSMISSION  COEFFICIENT  FOR  P  WAVES 


ooo  oo 
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...FROM  LAYER  1  TO  LAYER  2  IS  TPP12  =  PHID( 2 ) /PHID( 1 ) 

CX3MPLEX  GLl  ,GL2 
COMPLEX  RU.RD.TU.TO.CO 
COMPLEX  Y 

. . . NOTE  THAT :  GL= 1 /ALPHA 

DATA  C0/(0. 0,0.0)/ 

Y  s  RHOl*GL2  *  RH02*GL1 
TD  =  2.0*RHOl*GLl/Y 
RD  =  (RH02*GL1  -  RHOl*GL2)/Y 
TU  =  2.0*RHO2*GL2/Y 
RU  =  -RU 
RETURN 


SUBROUTINE  REVERB 


SUBROUTINE  REVERB(A.B) 

COOPLEX  A.B.X.Y.OET.CO.Cl 
LOGICAL  WARNED 

DATA  CO /(O. 0,0.0)/, Cl /(I. 0,0.0)/, WARNED/ . FALSE . / 

IF  (ABS(REAL(A) )«ABS(A1MAG(A) ).LT.  1.0E-12)THEN 
IF  (.NOT.  WARNED)THEN 

WRITE(6, ’ (/."•••MULTIPLES  IGNORED  WHERE  MAGNITUDE  i  E-12")’) 
WARNED  -  .TRUE. 

ENDIF 

...IGNORE  MULTIPLES,  TO  AVOID  UNDERFLOW  PROBLEMS 


B  =  Cl 
ELSE 

X  =  1.0  - 
Y  =  1.0 
DET  =  X^Y 
B  =  Y/DET 
ENDIF 
RETURN 


SUBROUTINE  FILTER 


I 


i 

r; 


V, 


■p.w. 
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SUBHOUTINE  FILTEK(REFW,NW.TSEC) 

LOGICAL  ^B>HASE 

COMMON  /FILTDr;  F](10),IDB(10),IFILTYP(10),NFILT,MPHASE, 

♦  FFl(lO) 


REAL  MAXABSV 

LOGICAL  FRSTTIM 

COMPLEX  FILT(1024),  FILT1(1024).  KEFW(1024) 

DATA  PI/3. 141592/ 

DATA  FRSTTIM/. TRUE./ 

IF  (NFILT  .EQ.  0)  RETURN 

C 

C  . GENERATE  FILTER . 

C 

IF  (FRSTTIM)  THEN 

C 

C  ...INITIATE  FILTER  IN  POSITIVE  FREQUENCY  DOMAIN.  NEGATIVE 

C  ...FREQUENCIES  BECOME  COMPLEX  CONJUGATES  OF  POSITIVE  FREQUENCIES 

C  ...AND  CAN  BE  COMPUTED  WHEN  NECESSARY. 

C 

DO  10  1=1,  NW+1 

FILT(I)  =  CMPLXd.O,  0.0) 

1 0  OONT I NUE 

C 

C  .. .SMEARING  RANGE 


C 

DELFREQ  =  0.0 

IF  (MPHASE)  DELFREQ=2*PI/TSEC 
DO  15  1=1,  NFILT 

IF  (IFILTYP(I)  .EQ.  1)  THEN 

C 

C  ...LOW  PASS  BUTTERWORTH  FILTER 

C 

POLES  =  FLOAT(IDB(l))/6  ♦  1 
R  *  l/((Fl(I)-DELFREQ)*TSEC) 

DO  20  J=l,  NW+1 

FILT(J)  =  FILT(J)/SQRT(l*{(J-l)*R)»»POLES) 
20  CONTINUE 

ELSEIF  (IFILTYP(I)  .EQ.  2)  THEN 

C 

C  ...HIGH  PASS  BUTTERWORTH  FILTER 

C 

POLES  =  FLOAT/ IOB(())/«  ♦  1 
FACTOR*  l/( (Fl( I )+DELFREQ)*TSEC) 

DO  25  J=l ,  NW+l 

R  =  ( (J-1 )*FACT0R)**P0LES 
FILT(4)  =  FILT(J)*SQRT(R/(1+R)) 

2  5  OONT I NUE 

ELSE 
C 

C  ...  NOTCH  F I LTER 

C 

IF( I .EQ.2 )GO  TO  26 

INDEXl  =  IFIX( (Fl( I )-DELFREQ)*TSEC)  +  1 
GO  TO  27 

26  INDEX2  =  IFIX((F1(I)*DELFR£Q)*TSEC)  ♦  2 


V 


'  V-:,-- '.s'.- 
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IF  (INDEX2-INDEX1  .LT.  16)  THEN 

V«11TE(6,'("  •••  Fl(l)  TOO  CLOSE  TO  Fl(2)  IN  NOTCH  - 

♦  "FILTERING  DONE  •••")•) 

RETURN 

ENUIF 

...IDB  REDUCTION 
00  TO  29 

27  R  -  1  -  10«*(-(FLQAT(IDB(I))/20)) 

. . .NEILS  DESIGN 
IFd.EQ.DOO  TO  28 

29  DO  30  J=IN0EX1,  INIKX2 
FILT(J)  =  FILT(J)* 

♦  (1  ♦  R*(OOS(2*PI*(J-INDEXl)/(INI»X2-INDEXl))-l)/2) 

30  CONTINUE  I 

28  ENDIF  ^ 

15  CONTINUE 

IF  (IVPHASE)  THEN 

. . .SET  A(W)  =  LN  F(W) 

DO  35  1=1,  NW^l 

IF  (REAL(FILT(I))  .EQ.  0.0)  THEN 
FILT(I)  -  CM>LX(-30.0,  0.0) 

ELSE 

FILT(I)  =  CM>LX(ALOG(REAL(FILT(I))),  0.0) 

ENDIF 

15  CONTINUE 

...  NEGATIVE  NEEDED  FC»  FFT 

DO  40  I«NW^2.  NW*2 

FILT(I)  »  OONJG(FILT(NW*2-I*2)) 

10  CONTINUE 

...A(T)  »  INVERSE  FFT  A(W) 

CALL  FFT(NW*2,FILT,-l.O) 

C  ...  BIT)  =  A(T)*G(T) 

DO  45  1=2,  NVKl 

C  ...  GIT) 

R  =  Il+OOSIPI*! I-l )/NW))  /2 
C  . . .  NW+ 1  VALUE  =  0 

FILTH)  =  R'FILTII) 

C  ...1  ST  VALUE  UNCHANGED 

FILTINW»2+2-I )  =  R*FILTINW*2»2-I) 

C  ...2  &  NW*2  VALUE  NEARLY  =  0 

45  CONTINUE 

C  ...NVKl  VAL  =  0  IS  USED  HERE 

DO  50  1=1,  NW 

FILTH  I)  »  FILTH) 

50  CONTINUE 
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C  ...  C(T)  HAS  INVERTED  2  HALF 

DO  55  l=NW*l,  NW*2 

FILTl(I)  =  -OONJG(FILT(l)) 

55  CONTINUE 

C  ...  B(W)  =  FT  B(T),  C(W)  =  FT  C(T) 

CALL  FFT(NW*2,FILT,1.0) 

CALL  FFT(NW»2,FILT1,1.0) 

MAXABSV  =0.0 

C  ...  F'(W)  =  EXP(B(W)*C(W)) 

DO  60  1=1,  NW^l 

FILT(I)=CEXP(FILT(I)  ♦  FILTl(I)) 

MAXABSV  =  AMAXl ( CABS (FILT( I)), MAXABSV) 

60  OONTINUE 

DO  65  1=1 ,  NW+1 

C  ...NORMALIZE 

FILT(I)  =  FILT( I )/MAXABSV 
65  OONTINUE 

END  IF 

FRSTTIM  =  .FALSE. 

ENDIF 


C 

C  . ACTUAL  FILTERING  FOLLOWS . 

C 

C  ...MULTIPLY  BY  POSITIVE  FREQ  FILTER 

C 


DO  70  1=1,  NW+1 

REFW( I )  =  REFW( I )  •  FILT(I) 

70  OONTINUE 
C 

C  ...  NEGATIVE  FREQUENCY  ARE  CONJUGATES 
C 

DO  75  I=NW+2,  NW*2 

REFW(I)  =  OONJG(REFW(NW*2-I*2) ) 

75  OONTINUE 
RETURN 
END 


0 . 

c 

O  SUBROUTINE  FFT 

0 

O . 

SUBROUTINE  FFT( LX,CX,S IGN I ) 
OOMPLEX  GX(LX) ,CTEMP,GW 
J  =  1 

SC  =  SQRT( 1 . /FLOAT! LX) ) 

DO  3  1=1 ,LX 
IF( I .GT. J )  GO  TO  1 
CTEMP  =  CX(J)*SC 
CX{J)  =  CX( I )*SC 
CX( I )  =  cn  EMP 

1  M  =  LX/2 

2  IF(J.LE.M)  GO  TO  3 


vv 


ooooo 
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J  s  J  -  M 
M  3  M/2 

IF(M.GE.l)  00  TO  2 

3  J  -  J  M 
L  =  1 

4  I  STEP  3  2*L 
DO  S  M»1,L 

AA  s  3141592SS353979.D-14*S1GNI«(M-1)/L 
CW  -  CM>LX(CX)S(AA).SIN(AA)) 

DO  5  1=M,LX, ISTEP 
CTEIi*  =  CW»CX(I+L) 

CX(I4^L)  >  CX(I)  -  CTEM* 

5  CX(I)  s  CX(I)  ♦  CTEM> 

I  s  IQTFP 

IF(L.LT.LX)  00  TO  4 

RETURN 

END 


SUBROUTINE  WRTPLT 


SUBROVTINE  WRTPLT 

OOKMON  /  NUM  /  CL( 200 ) ,QL( 200 ) .RHO( 200 ) .T( 200 ) ,EPS ,S ICMA, 

>  CRUST. NL.NTRACE. LOBS, NW,TSEC,TSC,ASC,NMULT, 

*  TLAG.YSIDE.TLNTH 
OOOMM  /TS/  REFT(1024) 

OdkMON /AR/ 1 AOC ,  WI NOOW 

OOMiON  /FILTDT/  Fl(  10  ) ,  IDB(  10  ) ,  IFILTYPI 10  )  .NFILT,M>HASE. 

♦  FFl(lO) 

COMMON  /  CHAR  /ITIME, lOATE, lOBT, TITLE 
CHARACTER* 1  lOBT 
CHARACTER*10  ITIME, IDATE 
CHARACrER*80  TITLE 

OPEN(UHIT»4,FILE»’PLTIN' ,STATUS=’NEW'  ) 

IFdOBT.EQ.  'P'  )THEN 
INUM^l 

ELSE  IFdOBT.EQ. 'V  )THEN 
INUM=2 
ELSE 
INUM>3 
END  IF 
NT*2*NW 

WRITE(4,1001)TITLE 
1001  FORMAT(A) 

WR I TE ( 4 , 9 9 9 ) NTRACE , YS I DE , NW , TLNTH . NL , NF 1 LT ,  M>HASE 
999  FC«MAT("NTRACE=", I4,",YSIDE»",F5.2,",NW«", 15, ".TLNTHs", 
♦F5.2,",HL*", I3,",NFILT»", I3,",**»HASE*",L2) 
WRITE(4,9I9)TSEC,TSC,TLA0,ASC 

989  FORMAT!  "TSEO  "  ,F10 . 0  ,  "  ,TSO "  ,F10 . 5 ,  " ,TLAG» "  ,F10 . 0  ,  " , ASC=  "  ,F10 . 5 ) 
DO  9  I^l.NFILT 


n n o  cooo  nonnn 


28 


GENERATION  OF  VERTICALLY  INCIDENT  SEISMOGRAMS 


WR1TE(4,99S)IFILTYP(I),I0B(I),FF1(I) 

995  FORMAT! "FILTER  TYPE* " , 12 , IX, "DB-SLOPE= " , 1 3 , IX, 

•"CUT-OFF  FREQUENCY* ",F5.1) 

5  CONTINUE 

WR1TE(4,987) INUM, LOBS, NMULT, CRUST, lAGC, WINDOW 
987  FORMAT! IX, "IOBT=", II , ",LOBS=", 13 , ",NMULT= " , 1 5 , " .CRUST* " , FlO . 5 
♦ , " , lAGC* " , 13 , " , WINDOW" ,F7 . 2 ) 

DO  6  1*1, NL 

WRITE! 4, ’ (IX, "CL=",F5.2,",QL=",F6.0,",RHO*",F5.2,",T»",F6.1)’ ) 
♦CL!I),QL!I),RHO!I),T!I) 

6  CONTINUE 

DO  10  1=1, NT 

WRITE!4, ' !"REFT=",1X,F20.14)* )REFT!I ) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  AGCFREQ 


SUBROUTINE  AGCFREQ 
DIMENSION  SQR£FT(1024) .BB!2) 

COMPLEX  SQREFW! 1024) .  WW!l024) 

COMMON  /  NUM  /  CL! 200 ) ,QL! 200 ) ,RHO( 200 ) ,T! 200 ) .EPS , S ICMA, 
CRUST , NL , NTRACE , LOBS , NW, TSEC , TSC , ASC , NMULT , 
TLAG.YSIDE.TLNTH 

COMMON  /  CHAR  /  ITIME, IDATE, lOBT, TITLE 
COMMON  /TS/  HEFT! 1024) 

OOMMON/AR/ lAGC, WINDOW 
DATA  PI/3. 141592653/ 

....SQUARE  TRACE  AMPLITUDES  AND  MAKE  COMPLEX 
....  INITIALIZE  FILTER  TO  COMPLEX  ZERO 

SUMSQ=0.0 
NT  *  NW»2 
DO  330  1  *  1,NT 

SQREFT!!)  =REFT(I)»REFT!I) 

SUMSQ=SUMSQ^SQREFT(  I ) 

SQREFW!!)  =CMPLX{ SQREFT! I ), 0 . 0 ) 

WW( I )  >  CMPLX!0.0,0.0) 

CONTINUE 

CONDITION  THE  INPUT  TIME  SERIES 

AVSQ=SUMSQ/NT 
BOOST* I . 0 /AVSQ 
DO  331  1=1, NT 

SQREFW!!)  =  SQREFW! i ) "BOOST 
331  CONTINUE 
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C 

C  . . .  .CXMSTRUCT  AOC  FILTER  IN  TIME  DGM4IN 

C  - "WINDOW"  (READ  IN  BY  USER)  IS  FILTER  LENGTH  IN  MSECS 

C 

DELTAT  »  TSEC/NT 

ITW2  «  NINT((WINDOW/DELTAT)/3.0) 

DO  340  I  >  1,ITW2>1 

m(l)  *CM>LX(0.5*(1.0«^CX)6(PI«(I-1)/1TW2)),0.0) 

J  ■  NT-^l-I 
WW(J)  >  WW(I) 

340  CONTINUE 

C 

C  CONDITION  THE  FILTER 

C 

SUMFIL^O.O 
DO  341  1=1, NT 

SIMF1L=SUMF1L  «  WW( I ) 

341  CONTINUE 

B06T  =  NT/SUMFIL 
DO  342  1=1, NT 

WW(  1 )=WW(  1 ) "BOOST 

342  CONTINUE 
C 

C  ...TRANSFORM  TO  FREQUENCY  DOMAIN  AND  FILTER 

C 

CALL  FFT(NT,SqREFW,^1.0) 

CALL  FFT(NT,WW,^1.0) 

DO  350  I  =1,NT 

S<^FW(I)  =  SQREFWI I  )«WW(  I ) 

350  CONTINUE 
C 

C  ...TRUNCATE  DC  OQM>ONENT  TO  5  SIGNIFICANT  FIGURES 
C 

TEM>=  REAL  ( SQR£FW(  1 ) ) 

ENCODE!  14 , 345  ,BB)TEMP 
raECXXS  ( 1 4 , 3  4  5 ,  BB )  TEM>  1 

345  F0R«MT(1E14.5) 

SQREFW(l)  =  CM>LX(TEM>  1,0.0) 

C 

C  ...RETURN  TO  TIME  DOMMN 
C 

CALL  FFT(NT,SQREFW, -1.0) 

C 

C  RESTORE  TIME  SERIES 

C 

SUMSQF  =0.0 
DO  346  1=1, NT 

SUMSQF  =  SUM5QF  ♦  S(^FW(I)**2 

346  CONTINUE 

AVSQF  =  SUMSQF /NT 
BOOST  =  AVSQ/AVSQF 
DO  347  1=1, NT 

SQREFW(I)  «  SQREFW( I )«BO06T 

347  CONTINUE 
C 


o o  o  o  o  no 
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APPLY  AGC 

CURRMIN  =  REAL(SQREFW{1)) 

DO  360  1=  2, NT 

CURRMIN  =  AMIN1(CURRM1N,R£AL(SQREFW(  I  ))) 

60  CONTINUE 

DO  380  I  =  I, NT 

IF  (CURRMIN. LT. 0.0)  THEN 

REFT(I)  =  REFT( I )/SQRT(REAL(SQREFW(  I ) )  +  ABS (CURRMIN*2 . 0 ) ) 
ELSE 

REFT(I)  =  REFT( I )/SQRT(REAL(SQREFW( I ) ) ) 

ENDIF 

80  CONTINUE 
RETURN 
END 


SUBROUTINE  TERM  IN 


SUBROUTINE  TERMIN(NERR,  LNR) 


COMVION  /  NUM  /  CL(200  )  ,QL(200)  ,RHO(200  )  ,T(200  )  , EPS, SIGMA, 

♦  CRUST , N  L , NTRACE , LOBS , NW , TS  EC , TSC , ASC , NMULT , 

♦  TLAG,YS IDE, TENTH 
COMMON  /  CHAR  /  ITIME, IDATE, lOBT, TITLE 

IF  (NERR  .LT.  10)  THEN 

WRITE(6,'(/"  •••  ERROR  IN  MODEL  FILE  LAYER:  ",I4)')  LNR 
ELSEIF  (NERR  .LT.  20)  THEN 
CALL  PRTMLDA 

WRITE(6,'(/"  •••  ERROR  IN  PIUXESSED  LAYER:  ",I4)')  LNR 
ENDIF 

IF  (NERR  .EQ.  I)  THEN 

WRITE(6,’("  TOO  MANY  LAYERS  IN  MODEL:  ",I6,”,  MAX=200")’)  NL 
ELSEIF  (NERR  . EQ.  3)  THEN 

WRITE(6, ' ( "CL  CANNOT  BE  0  AND  H  CANNOT  BE  -VE  FOR  FIRST 
♦LAYER" ) ' ) 

ELSEIF  (NERR  .EQ.  4)  THEN 

WRITE(6, ' ( "LAST  LAYER  CANNOT  BE  AN  INTERPOLATED  LAYER")') 

ELSEIF  (NERR  . EQ.  S)  THEN 

WR I TE( 6, '( "ABSOLUTE  DEPTH  GIVEN  IS  LESS  THAN  THE  CURRENT  DEPTH")') 

ELSEIF  (NERR  .EQ.  1 1 )THEN 

WRITE(6, ' ( "CL  NOT  IN  RANGE  0  i  CL  i  10")') 

ELSEIF  (NERR  . EQ.  I2)THEN 
WRITE(6,  '("QL  i  100")  '  ) 

ELSEIF  (NERR  . EQ.  15)  THEN 

WRITE(6, ’ ("RHO  NOT  IN  RANGE  0  i  RHO  i  10")') 

ELSEIF  (NERR  . EQ.  21)  THEN 

WHITE(6,'(/"  •••  TOO  FEW  (OR  TOO  MANY)  VALUES  ON  ONE  LINE  " 

♦  "OF  MODEL  FILE  CHECK  FOR  /  AT  END  OF  INTERPOLATED  LAYER  ", 
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Test  3.  Multiples 


Two  layers,  75  and  100  m  thick,  are 
sandwiched  between  half  spaces  (fig.  3B).  Q  = 
200  and  density  =  1.1  for  all  units.  The 
100-meter-thick  layer  has  a  velocity  of  3.0 
km/sec;  the  other  units  have  a  velocity  of  1.5 
km/sec.  The  seismogram  consisted  of  the 
time-domain  reflection  coefficient  plotted 
without  filtering. 

A  series  of  tests  were  run  with  options 
beginning  with  “primaries  only”  to  options 
for  a  selected  number  of  “multiples”  and 
finally  a  “complete-response”  seismogram 
(fig.  5).  R1  and  R2  identify  the  reflections 
from  the  top  of  the  two  reflectors.  Ml 


through  M4  identify  the  “peg-leg”  'lUll^ples. 
Figure  5A  demonstrates  the  ..omplete- 
response  option;  the  clipping  occurs  because 
of  a  large  amplification  option  (ASC  =  100). 
Figure  5B  shows  the  seismogram  for  the  same 
model  but  for  a  user  option  of  only  one 
multiple.  Note  that  no  reflections  are 
observed  beyond  the  first  multiple  arrival  at 
234  ms.  Arrival  times,  amplitudes,  and 
polarity  in  figure  5  are  consistent  with  those 
predicted  by  hand  calculations  based  on 
reflection-transmission  coefficients.  Table  3 
summarizes  the  results. 


Table  3.  Arrival  times  of  reflections  and  multiples  for 
a  two-layer  model  (fig.  5) 


Reflector 

index 

Primary 

only 

One 

multiple 

Two 

multiples 

Full 

response 

R1 

100  ms 

100  ms 

100  ms 

100  ms 

R2 

167  ms 

167  ms 

167  ms 

167  ms 

Ml 

- 

234  ms 

234  ms 

234  ms 

M2 

- 

— 

300  ms 

300  ms 

M3 

— 

— 

- 

367  ms 

M4 

— 

— 

— 

434  ms 

Table  4.  Input  data  used  to  test  a  full  response  (fig.  5A) 


TEST 

3  -  D2L. 

FULL  RESPONSE  OPTION. 

EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6. 

CL 

QL 

RHO 

THICKNESS 

1.50 

200. 

1.1 

0. 

1.50 

200  . 

1.1 

75. 

3.00 

200. 

1.1 

100. 

1.5 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC 

#FREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL 

AGC 

WINDOW 

5 

512 

1000. 

006  100.0  -1 

0 

0 

DETECTOR  TYPE 

DETECTOR 

LOCATION 

REFLECTION  LAYER  1 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER-TYPE  OB  SLOPE  CUT-OFF  FREQUENCY 
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TEST  2  -  OIL.  HIGH  Q  (=1500) 

reflection  coefficient 

AT  0B5  IN  LAYER  1  AT  .0  H  DEPTH 
FREQtNWl=  512  T ( MSECS  1= 2000 

MIN  =  .50  HZ  MAX  =256.00  HZ 

PRIMARIES  ONLY 
NO  A5C.  ASC  =  1. 


TEST  2  -  DIL.  LOW  Q  (=5) 

REFLECTION  COEFFICIENT 

AT  OB5  IN  LAYER  1  AT  .0  M  DEPTH 

FREQ(NW)=  512  T ( MSECS )= 2000 

MIN  =  .50  HZ  MAX  =256.00  HZ 

PRIMARIES  ONLY 

NO  AGC.  ASC  =  1. 


NO  filters 


NO  FILTERS 


GENERATION  OF  VERTICALLY  INCIDENT  SEISMOGRAMS 
Test  2.  Attenuation  Factor  Q 


A  layer  1000  m  thick  was  overlain  by  a  half 
space  with  an  identical  velocity  (1.5  km/sec) 
and  underlain  by  a  half  space  with  a  velocity 
of  3.0  km/sec  (fig.  3A).  All  units  had  a 
constant  density  of  1.1  gm/cc.  No  filter  was 
applied.  The  seismograms  consisted  of  a  plot 
of  reflection  coefficients  for  primary  arrivals 
only.  A  wide  frequency  range  of  1-512  hz  was 
invoked  to  study  the  waveforms  better. 

The  Q  of  rocks  for  seismic  energy  is  of  the 
order  50,  althou^  poorly  consolidated 
materials  may  be  lower.  In  this  test  Q  was 
varied  from  5  to  2000.  For  Q  values  above  50, 
the  response  was  essentially  a  spike  with  a 
width  of  less  than  10  msec.  Figure  4A 
illustrates  a  Q  of  1500.  For  Q  values  of  20  or 
less,  the  pulse  width  of  the  response  increased 


to  20  msec,  with  a  width  of  about  150  msec 
for  a  Q  of  5  (fig.  4B).  As  expected,  low  Q 
emphasized  lower  frequencies. 

A  change  of  Q  from  15  in  the  top  half 
space  to  1500  in  a  layer  1,000  m  thick  (but 
with  no  change  in  the  user-specified  velocity 
between  these  two  regions)  yielded  a  small 
pulse  at  zero  time.  (See  arrow  in  fig.  4A.)  This 
is  consistent  with  our  algorithm,  which  brings 
in  attenuation  changes  through  a  complex 
velocity. 

The  following  records  (table  2)  were  used 
as  input  to  Program  VISP  to  generate  figure 
4A.  The  same  records  were  used  to  generate 
figure  4B,  except  that  QL  values  were  all  set 
equal  to  five. 


Table  2.  Input  data  used  to  test  high  Q  (fig.  4A) 

TEST  2  -  DIL.  HIGH  Q  (=1500) 

EPS  SIGMA  YSIDE  XSIDE 

0.001  0.1  40.  6.0 

CL  QL  RHO  THICKNESS 

1.5  15.  1.1  0. 

1.5  1500.  1.1  1000. 

3.0  1500.  1.1  0. 

9999  9999  9999  9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  2000.  .003  1.0  000 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER-TYPE  DB  SLOPE  CITT-OFF  FREQUENCY 
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Appendix  4.  Input  Records  and  Output  Plots  for  12  Tests 

Thirteen  tests  were  run  to  demonstrate  the  following  test  cases  is  related  to  one  of  these 

reliability  of  the  options  available  to  the  user  three  models,  although  the  velocity,  Q,  and 

of  Program  VISP.  Test  1  consisted  of  a  model  (or)  density  may  be  varied.  Program  VISP  has 

simulating  a  standard  field  record  (fig.  2).  In  detailed  comment  cards  at  the  beginning  that 

the  remaining  12  tests,  three  basic  models  explain  the  various  user  options  tested  in  this 

were  used:  (1)  one-layer,  (2)  two-layer,  and  appendix. 

(3)  six-layer  (figs.  3A-3C).  Each  of  the 


CL  IKM/SECl  CL  I  KM/SEC)  CL  tKfl/SEC) 


1  2  3  1  2  3  1  2  3 


1  2  3  1  2  3  1  2  3 

Figure  3.  Velocity  models  for  TesU  2-13  (appendix  4).  A,  one-Uyer;  B,  two-layer; 

C,  six-layer. 


noon 


GENERATION  OF  VERTICALLY  INCIDENT  SEISMOGRAMS 
IF  (MODdXRANGE.lOO)  .EQ.  0) 

+  CALL  NUNBER(X-0.15,Y+YOFFSET,0.055,FLOAT(IXRANGE),0.0,-1) 

CALL  PLOT(X,Y, INKOFF) 

CALL  PLOT(X+0.75,Y, INKON) 

IF  (I.EQ.O)  CALL  SYNBOL(X0>TLNTH/2-0 . 4 ,  Y-0 . 7 , 0 . 16  , 

+  "METERS ”,0.0,6) 

Y  =  Y  +  7.0  -  0.4 
X  =  XO 

PLOT  MODEL  PARAMETERS 


DO  170  1=1,3 

IF  (PMAXd)  .NE.  PMINd))QO  TO  55 
YSC  =  0.0 
Y1  =  YSTRIP/2 
GO  TO  65 

55  YSC  =  YSTRIP/(PMAXd)-PMINd)) 

Y1  =0.0 

65  YI  =  YO  ♦  0.2  +  d-l)*(ySTRIP+0.2)  +  0.1 
IF  (d.EQ.  1)  .OR.  d.EQ.5)) 

♦  CALL  PLOT(X0,YI+(PARMd,I)-PMINd))*YSC+Yl,  INKOFF) 

X  =  XO+0.75 

LIMIT=NL-1 

DO  180  J=l, LIMIT 

IF  ((J.NE.l)  .OR.  (d.LT.2)  .OR.  (I  .GT.  4  ) )  )GO  TO  75 
CALL  PLOT(X,YI+(PARM(2, I)-PMIN( I ) )*YSC+Y1 , INKOFF) 

GO  TO  85 

75  CALL  PLOT{X,YI  +  (PARM(J,I)-PMINd))*YSC+Yl,INKON) 

CALL  PLOT(X,YI  +  (PARM(J  +  l,I)-PMINd))*YSC+Yl,  INKON) 

85  IF  (J.LT.NL-1)  X=X  +  T( J+1 )*(TLNTH-1 . 5 ) / IXRANGE 

180  CONTINUE 

X  =  XO+TLNTH 

CALL  PLOT(X,YI+(PARM(NL,  I)-PMINd)>*ySC+yi,  INKON) 

I F (I . EQ.  1 ) CALL  SYMBOL ( XO -0 . 6 . YO  +  d - 1 ) • ( YSTR IP+.2)  +  .3,.l, 
+  "CL  (KM/SEC)  ”,90., 12) 

IFd.EQ.2)CALL  SY^BOL(  XO -0 . 6  ,  Y0+ (I -1 )  •  ( YSTR1P+ .  2  )+ .  3  , .  1 , 
+"QL*  ”,90., 12) 

IFd  .EQ.3)CALL  SYNBOL(X0 -0 . 6 ,  Y0  + d -1 )  *  ( YSTRIP+ .  2  )  +  .  3  , .  1 , 
+"RHO  (GM/OC)  ”,90. ,12) 

CALL  NUNBER( X0-0.6,YO  +  I.96,0.10, PARMP Y, 90.0,3) 

170  CONTINUE 

XO  =  XO+TLNTH+2.5 

RETURN 

END 
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10  CALL  PLOT(X,Y.  INKOFF) 

IF(J.EQ.2)CALL  PLOT(X+.  !•(  1*2-1 )  ,Y,  INKW) 

CALL  PLOT(X+O.l*(l*2-1),Y,INKON) 

CALL  NUAeER(X-^»»'FSET,Y.0.1,FLQAT(PMIN(J)),90.0,-l 
CALL  PLOT! X,Y, INKOFF) 

ISTART=PMIN(J)>1 

ISTOP=PMAX(J) 

00  130  K=  I START, I STOP 
DO  140  Lsl.lO 

Y  =  Y  +  YSTRIP/((PMAX(J)-PMIN(J))*10) 

CALL  PLOT(X,Y, INKON) 

IF(J.NE.2)  CALL  PLOT(X+0.05*( 1*2-1 ) ,Y, INKON) 
CALL  PLOT{X,Y, INKOFF) 

140  CXOTINUE 

CALL  PLOT(X+O.l*(  1*2-1  ),Y,INK<») 

IF(J.NE.2) 

+  CALL  NUIVBER(X+XOFFSET,Y-0.1,0.1,FLOAT(K),90 

IF(J.EQ.2.AND.K.EQ.PMAX(2)) 

♦  CALL  NUMBER(X+XOFFSET,Y-0.1,0.1,FLQAT(K),90 

CALL  PLOT(X.Y, INKOFF) 

130  CONTINUE 

20  Y  =  Y+0.1 

120  CONTINUE 

X  =  X+TLNTH 
Y  =  YO  +  0.2 
no  CONTINUE 
C 

C  PLOT  X  AXIS 

C  . 

C 

X  =  XO 

Y  s  YO  +  0.2 
I  *  0 

YOFFSET  =  -0.1+0.25*(I*2-1) 

CALL  PLOTCX.Y, INKOFF) 

X  =  X+0.75 

CALL  PLOT(X,Y, INKON) 

DO  160  N=1,IXRANGE,10 
J=N-1 

IF(MOD(J, 100). NE. 0)00  TO  30 

CALL  PLOT(X,Y>0.1*( 1*2-1), INKON) 

CALL  NUNBER( X-0. 05, Y* YOFFSET, 0. 055, FLOAT(J),0. 0,-1 
00  TO  40 

30  CALL  PLC7r(X,Y*0.05*(I*2-l), INKON) 

40  CALL  PLar(X,Y, INKOFF) 

X  =  X  +  (TLNTH-1.5)/IXRANGE*10 

CALL  PLOT  (X,Y, INKON) 

160  CONTINUE 

CALL  PLOT(X,Y+0.1*( 1*2-1 ),INKCM<) 
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INKON=2 

INKOFF=3 

Y0=0.5 

XO  =  XO+2.0 

PARMAX  =  ABSPARM(1,2) 

DO  50  1=  2,NL 

PARMAX  =  AMAX1(PARMAX,ABSPARM(I,2)) 

50  CONTINUE 

1F(PARMAX.LT. 100.0)  PARMPY  =  l.O 

IF(PARMAX.GE. 100.0. AND. PARMAX.lt. 1000.0)  PARMPY  =  0.1 
IF(PARMAX.GE. 1000.0. AND. PARMAX.lt. 10000.0)  PARMPY  =  0.01 
IF(PARMAX.GE. 10000.0)  PARMPY  =  0.001 
DO  60  1=1, NL 
DO  70  J=l,3,2 

PARM(I.J)  =  ABSPARMd.J) 

70  CONTINUE 

PARM(I,2)  =  ABSPARMl I ,2)*PARMPY 
60  CONTINUE 

DO  80  1=1,3 
PMAXd)  =  0 
PMlN(l)  =  1000 
80  CONTINUE 

PMlN(l)  =  PARM(1,1) 

DO  90  1=2, NL 
DO  100  J=l,3 

PMAX(J)  =  MAXl(FLOAT(PMAX(J)),PARM( I,J)+0.999) 
PMIN(J)  =  MINI ( FLOAT (PMIN(J)),PARM(  I, J)) 

100  CONTINUE 
90  CONTINUE 

YSTRIP  =  1.0 
IXRANGE  =  CRUST+0.99 

PLOT  Y  AXIS 


X  =  XO 

Y  =  YO  +  0.2 
DO  110  M=l,2 
I=M-1 

XOFFSET  =  0. 12  +  0. 3*(  1*2-1) 

CALL  PLOT(X, Y, INKOFF) 

DO  120  J=l,3 
Y  =  Y+0.1 

IF(PM1N(J).NE.PMAX(J))00  TO  10 

Y  =  Y  +  YSTRIP/2 
CALL  PLOT! X,Y, INKOFF) 

CALL  PLOT(X+0. 1*(  1*2-1  ),Y,  INKW) 

CALL  NUMBER! X+ XOFFSET, Y, 0.1, FLOAT(PMIN(J)), 90. 0,-1) 

Y  =  Y  +  YSTRIP/2 
GO  TO  20 
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C  PLOT  TIME  SERIES 

C  . 

c 

CURRMAX  »  ABS(REFT(1)) 

NT=NW*2 

DO  290  1=2, NT 

CURRMAX  =  AMAX1(CURRMAX.ABS(REFT(I))) 

290  CX)NTINUE 
X  =  XO 

y  =  YO  -  2.0 

c 

c  EFFECTIVELY  A  NULL  TRACE  (NOISE) 

C 

IF  (CURRMAX.LT.0.1)  CURRMAX=-10.0 
DO  300  ITHETA=1,NTRACE 
AY  =  (ITHETA-0.5)*XSC 
CALL  PLOT(X,Y*AY, INKOFF) 

C 

C  ALLOWS  FOR  TRUNCATED  TIME  SERIES  (TSC  i  0) 

C 

NPTS  =  NW*2*AMIN1(1.0,YSIDE/(TSEC*TSC)) 

YSC  =  NPTS*TSEC*TSC/( (NPTS-1)*NW*2) 

DO  310  1=1, NPTS 

1F(.N0T.  ASC.EQ.1.0)  00  TO  S5 
XX(I)  =  X+(I-1)*YSC 
YY(I)  =  Y+AY+REFT(I)/(CURRMAX*2) 

CALL  PLOT(XX(I),yY(I),INK(W) 

00  TO  310 

55  XI  =  AM1N1(ABS(REFT(I))*ASC/(CURRMAX*2),  0.4) 

IF  (REFT(I)  .LT.  0.0)  X1=-X1 
XX(I)  =  X+(I-1)*YSC 
YY(I)  =  Y+AY+Xl 
CALL  PLOT(XX(I),YY(I),INKCW) 

310  CONTINUE 
300  CONTINUE 
RETURN 
END 


SUBROUTINE  PARMPLT 


SUBROUTINE  PARMPLT 

OOMWON/  PLTPAR  /NTRACE , YS IDE, TSEC,TSC,TLAG,ASC,REFT( 1024 ) ,NW 
♦, lOBT, LOBS, NMULT,NL,CL( 200), ^(200), RHO( 200), T( 200), TLNTH,X0 
♦.CRUST, lAGC, WINDOW, NFILT 

CCXVMON/  FILTDT  /  IDB(  1 0  ) ,  IF ILTYP(  10 )  ,*i»HASE, FFl  ( 10 ) 

REAL  ABSPARM(200,3),  PARM(200,3),  TLNTH 

EQUIVALENCE  (ABSPARM( 1 , 1 ) ,CL( 1 ) ) 

1 NTEGE.'  PMAX(  3  ) ,  PMI N ( 3  ) 
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C 

X  =  XO 

Y  =  YO  -  2.4 
1  =  0 

XOFFSET  =  0. 07+0. 3*( 1*2-1) 

CALL  PLOT(X,Y, INKOFF) 

DO  270  M=1,IYRANGE 
J=M-1 

IF(.NOT.MOD(J-l,MODF)  .EQ.  INT(TLAG)-l)  GO  TO  50 
CALL  PLOT(X,Y+0.1*(I*2-1),INK(») 

C 

C  CHANGE  SIZE  OF  CHARACTER 
C 

IF(TSC.LT. ( .01 ) )00  TO  49 

CALL  NUMBER(X-0.07,Y+XOFFSET,0.055,FLQAT( J)-TLAG,0. , -1 ) 

GO  TO  50 

49  CALL  Ni;iVBER(X+.03,(Y+XOFFSET)-.05,0.055,FLOAT(J)-TLAG,90.  ,-l) 

50  1F(.NOT.MOD(J,(MODF-1)/10+1).EQ.O)  GO  TO  51 

CALL  PLOT(X,Y+0.07*( 1*2-1), INKOI) 

51  CALL  PLOT(X,Y, INKOFF) 

IF( .NOT. lYRANGE.LE. 20)00  TO  52 
DO  280  K=l,10 
Y  s  y  ♦  T^P*0  1 

CALL  PLOTIX.yI INKON) 

CALL  PLOT(X,Y+0.05*{ 1*2-1), INKON) 

CALL  PLOT! X.Y, INKOFF) 

280  CONTINUE 

GO  TO  270 

Y  =  Y+TSP 

CALL  PLOT(X,Y, INKON) 

270  CONTINUE 

CALL  PLOT(X,Y+0.1*(I*2-1),INKC»J) 

IFITSC.LT. ( .01) )00  TO  53 

CALL  NU5«ER( X-0 . 0 7 , Y+  XOFFSET , . 0 5  5 , FLOAT ( I YRANGE ) -TLAG , 0 . 0 , - 1 ) 

GO  TO  54 

53  CALL  NUMBER(X+.022,(Y+XOFFSET)-.05,.055,FLOAT(IYRANGE)-TLAG, 

*90. ,-l) 

54  Y  =  YO  -  1.4 

X  =  X0+IYRANGE*TSC/2-1.5 

CALL  SYMB0L(X+.75,Y-1.5,.15,"TIME  IN  MSEC  ",0.0,13) 

C 

C  PLOT  X  AXIS 

C  . 

C 

Y  =  YO  -  2.0 

CALL  PLOT! X0,Y, INKOFF) 

CALL  PLOT(X0,Y+0.4+YS, INKON) 

C 
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33  CALL  SY^©OL(XO♦1.8,6.5..1,"— filters  =  ZERO  PHASE---\90.  ,26) 

34  DO  101  I=1,NFILT 
lF(IFILTYP(I).NE.l)aO  TO  35 

CALL  NI»BER(X0+1.9+I*.15,7.5,.1,FF1(I),90. ,0) 

CALL  NUMBER(X0+1.9+I*.15,8.3,.1,FLQAT(IDB(I)),90. ,-l) 

CALL  SY5®OL(XO+1.9^I*.15,6.5,.l, 

•"HIGH  COT  AT  HZ;  DB  SLOPE" ,90 33 ) 

35  IF(IFILTYP(I).NE.2)aO  TO  36 

CALL  N05«ER(XO+1.9+1*.15,7.5,.1,FF1(1),^.  ,0) 

CALL  NOI«ER(X0+1.9+I*.15,8.3,.l,FLQAT(IDB(I)),90. ,-l) 

CALL  SY5«OL(XO+1. 9+1*. 15,6.5,.!, 

•"LOW  COT  AT  HZ;  DB  SLOPE" , 90 . , 33 ) 

36  IF({IF1LTYP(I).NE.1).AND. (IFILTYP(I).NE.2))00  TO  37 
GO  TO  101 

37  IF(l.EQ.2)GO  TO  101 

CALL  SYMBOL(XO+1.9+I*.15,6.5,.1, 

•"NOTCH  DB  FREQUENCY  RANGE  -  HZ", 90., 95) 


CALL  NU5BER(X0+1.9+I*.15,8.9,.l,FFl(l),9O. , 
CALL  NUMBER(X0+1.9+I*.15,9.4,.l,FFl{2),90. , 
CONTINUE 

IF(  lAGC.NE.DOO  TO  12 

CALL  SY\BOL(XO+1.4,6.5,.1,"AGC  WINDOW  * 
CALL  NUMBER(X0+1.4,7.8, .1 ,WINDOW,90. ,0) 
GO  TO  30 

CALL  SY5BOL(XO+1.4,6.5, .l,"NO  AGC.  ASC  = 
CALL  NUNBER(X0+1.4,7.9, .1 ,ASC,90. ,0) 

XO  =  XO+3.0 
XSC  =0.30 

FOR  TRUNCATED  TIME  SERIES 

lYRANGE  s  AMIN1(TSEC,YSIDE/TSC)  -  0.01  +  1 


MSECS", 90. ,24) 


',90. ,14) 


(TSC  i  0) 

IF  (lYRANGE  .LE.  10)  MDDF=1 
IF  ( lYRANGE. LE. 20 )MODF=2 
IF( lYRANGE. LE.50)MODF=5 
IF( IYRANGE.LE.100)MODF=10 
IFdYRANGE.CTT.lOO)  MODF=50 
YS  =  NTRACE"XSC 

PLOT  Y(TIME)  AXES 
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OC»MON/  FILTDT  / IDB( 10 ) , IFILTYP( 10 ) .MPHASE, FFl ( 10 ) 

DIMENSION  XX(I024),YY(1024) 

LOGICAL  MPHASE 

INKCXI=2 

INKOFF=3 

R=0.0 

Y0=9.0 

X0=0.75 

IF(  lOBT.EQ.  DCALL  SYNBOL(XO,6.S,  .1 , 

♦"PRESSURE  ",90.,24) 

IF( IOBT.EQ.2)CALL  SY\BOL( XO ,6 . 5 , . I , 

+  "VERTICAL  DISPLACEMENT  ",90.,24) 

IFdOBT.NE.l.AND.  lOBT.NE.  2  )CALL  SYACOLI XO ,  6 . 5 , 

+.1, "REFLECTION  COEFFICIENT  ",90.,22) 

IF(L(»S.NE.O)QO  TO  10 
GO  TO  15 

10  CALL  SY5BOL(XO+,26,6.5,.l, 

♦"AT  OBS  IN  LAYER  AT  M  DEPTH  ",90.,39) 

CALL  NU5eER(X0^.26,8.1, .l.FLOAT(LOBS),90. ,-l) 

CALL  SY\BOL(X0^.52,6.5,.1,"FREQ(NW)=",90. ,9) 

CALL  NU5BER(X0^.52,7.4,.1,NW,90.,-1) 

CALL  SY5fBOL(X0^.52,8.5,  .1,"T(MSECS)  =  ",90.  ,9) 

CALL  Nl»BER(X0^.52,9.3, .1,TSEC,90. ,-l) 

K=LOBS-l 

IF<K.LT.2)  GO  TO  11 
DO  20  1=2, K 
R=R^T( I) 

20  CONTINUE 

11  CALL  NUNBER(X0^.28,8.7,.1,R,90.,1) 

15  CALL  SY^BOL(X0♦. 8,6.4, .1,"  MIN  =  HZ  ",90.,17) 

CALL  NUM3ER(X0^.8,7.2, .l,(l/TSEC)*1000. ,90. ,2) 

CALL  SYNBOL(X0^.8,8.l, .1,"MAX  s  HZ", 90. ,14) 

CALL  NUNBER(X0^.8,8.55, .l,(NW/TSEC)*1000. ,90. ,2) 
IF{NMULT.LT.0)CALL  SY5O0L( XO^I . 1 , 6 . 5 , . 1 , 

♦"FULL  RESPONSE  (ALL  MULTIPLES)  ",90., 29) 

IF(NMULT.EQ.0)CALL  SYMBOL( XO^l . 1 , 6 . 5 ,. 1 , "PRIMARIES  ONLY  ", 

♦90. ,14) 

IF(NMULT.GT.0)GO  TO  25 
GO  TO  28 

25  CALL  SYMBOL(X0^1,l,6.5,.l,"NU5eER  OF  MULTIPLES  =  ",90.,21) 

CALL  NUMBER(X0^1.1,8.7,0.l,FLOAT(NMULT),90. ,-l) 

28  IF(NFILT.GT.0)GO  TO  27 

CALL  SYMBOL! 2. 5, 6. 5,. 1, "NO  FILTERS" , 90 ., 10 ) 

GO  TO  31 

27  IF( .NOT.MPHASE)GO  TO  33 

CALL  SYMBOL(X0^1.8,6.5,.1,"---FILTERS  =  MINIMUM  PHASE---" , 90 . , 29 ) 
GO  TO  34 


ooooo  ooooo 
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989  FORMAT(SX,F10.0,SX,P10.5,6X,F10.0,SX,F10.5) 

IF(NFILT.EQ.O)QO  TO  985 
DO  15  1-1,NFILT 

READd  ,986 )  IFILTY?(  I  )  ,  IDB(  I )  .FFl  (  I ) 

986  FC«MAT(12X, 12,10X, 13,19X,F5.1) 

15  CONTINUE 

985  READd, 987)  I(»T, LOBS, NMULT, CRUST,  IAOC,WlNDOIV 
987  FORMAT(6X, d,6X, I3,7X, IS,7X,F10.5,6X, I3.7X,F7.2) 
NT=2*NW 
DO  5  1=1, NL 

READd  ,  102  )CL(  I )  ,(^(  I)  ,RHO(  I )  ,T(  I ) 

102  FORMAT ( 4X, F5 . 2 , 4X, F6 . 0 , 5X, FS . 2 , 3X, F6 . 1 ) 

5  CONTINUE 

DO  10  1=1, NT 

READd,101)REFT(I) 

101  FORMAT(5X,F20.14) 

10  CC84TINUE 
RETURN 
END 


SUBROUTINE  INITPLT 


SUBROUTINE  INITPLT 
COMMON/  CHAR  /TITLE(8) 

COMMON/  PLTPAR  /NTRACE,YSIDE,TSEC,TSC,TLAG,ASC,REFTd024  )  ,NW 
lOBT, LOBS, NMULT, NL,CL( 200), QL( 200), RHO( 200), T( 200), TLNTH,X0 
+, CRUST, IAGC,WINDOW,NFILT 

CO5M0N/  FILTDT  /  IDBdO  ) ,  IFILTYPdO  )  ,M>HASE, FFl  (10 ) 

CALL  SY5eOL(.S,.5,0.15,"VISP;  ",90. ,6) 

CALL  SY5BOL(.5,,8,.15,"  PLOTS  OF  PARAMETERS  ", 

^90.0,34) 

CALL  SYMBOLd.S,.5,.I,T(TLE,90.,80) 

CALL  SYAeOL(.5,6.5,.12,TITLE,90. ,80) 

X0=1.5 

RETURN 

END 


SUBROUTINE  RSPMDL 


SUBROUTINE  RSPMn. 

COMMON/  PLTPAR  /NTRACE,  YS  IDE,TSEC,TSC,TLAG,ASC,REFTd024  )  ,NW 
♦, lOBT, LOBS, NMULT, NL,CL( 200), ^(200), RHO(  200), T( 200), TLNTH,X0 
♦ , CRUST , I AGC ,  WI NDOW, NF I LT 
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C 

C  PROGRAM  PLTVISP 

C 

C  (PLOTTING  ROUTINE  FOR  PROGRAM  VISP) 

C 

C  MODIFIED  FOR  CDC  855  USE  AT  INDIANA  UNIVERSITY 

C  FROM  PLOTTING  CODE  DEVELOPED  BY  FRAZER, BATES 

C  AND  RUDMAN  AT  THE  HAWAII  INSTITUTE  OF  GEOPHYSICS 

C 

C  INPUT  FOR  THIS  PROGRAM  IS  GENERATED  AS  A  FILE 

C  ON  TAPE  4  (AND  DESIGNATED  "PLTIN")  BY  PROGRAM 

C  VISP  (SEE  SUBROUTINE  WRTPLT). 

C 


MAIN  PROGRAM 


PROGRAM  TSTPLT(TAPE1, OUTPUT, TAPE4=OUTPUT, PLOT, TAPE7=PLOT) 
COMMON/  CHAR  /TITLE(8) 

COMMOI/  PLTPAR  /NTRACE ,YS IDE, TSEC,TSC,TLAG, ASC, REFT( 1024 ) ,NW 
lOBT, LOBS, NMULT,NL,CL( 200), QL( 200), RHO( 200), T( 200) ,TLNTH,X0 
♦, CRUST, IAGC,WlNDOW,NFILT 

COMVDN/  FILTDT  / IDB( 10 ) , IFILTYP( 10 ) ,MPHASE ,FFl ( 10 ) 

CALL  IDENT(7) 

CALL  READINP 

CALL  INITPLT 

CALL  PARMPLT 

CALL  RSPMDL 

CALL  PLOT(0. ,0. ,999) 

STOP 

END 

C . 

c 

C  SUBROUTINE  READINP 

C 

C . 

SUBROUTINE  READINP 
COMMON/  CHAR  /TITLE(8) 

COftMON/  PLTPAR  /NTRACE, YS IDE,TSEC, TSC,TLAG, ASC, REFT( 1 024 ) ,NW 
♦  , lOBT, LOBS, NMULT,NL,CL( 200),  <^(200), RHO( 200) ,T( 200) ,TLNTH,X0 
♦, CRUST, lAGC, WINDOW, NF I LT 

COMVDN/  FILTDT  / IDB( 10 ) , IFILTYP( 10 ) ,MPHASE , FFl ( 10 ) 
READ(1,1001)TITLE 
1001  FORMAT(8A10) 

READ( 1 , 9 9 9 ) NTRACE , YS I DE , NW, TLNTH , NL , NF I LT , MPHASE 
999  FORMAT! 7  X, 1 4 , 7X, F5 . 2 , 4X, 1 5 , 7X, F5 . 2 , 4X, 1 3 , 7X, 1 3 , 8X, L2 ) 

READ (  1 , 9  8  9  ) TSEC , TSC , TLAG , ASC 
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*  "DECLARATION  LINES")') 

ELSEIF  (NERR  .EQ.  23)  THEN  «  i*  ■  uiteT  RF  A  POMIER  OF  2" 

VfllITE(6.’(/"  •••  ERR<»  -  NFREQ  ",U.  »WST  BE  A  POWER  ur  z 

*  ",  1«-512")‘)  LNR 


)')  LNR 


ELSEIF  (NERR  .EQ.  30)  THEN  BprvTFRTFn.  « 

WRITE(6  '(/"  ERROR  -  1  lO  FILTERS  REQUESTED.  ) 

*  Lm 

ELSEIF  (NERR  .EQ.  31)  THEN  .  i  ,  i  MiifcKFR" 

WRITE(6  '(/"  •••  Layer  fwi  obs  placement  I  2  or  t  nu*©er  , 

♦  ’  "OF  LAYERS")' ) 

^^^1TE(6, '  C  •••  UNKNOWN  ERROR,  NU^«ER:  ",14)')  NERR 

ENDIF 

STOP 

END 
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ELSEIF  (NERR 

WRITE(6, ' (/ 

* 

+■ 

ELSEIF  (NERR 

WR1TE(6, '(/ 

■1 

ELSEIF  (NERR 

■ 

WRITE(6. '(/ 

hi' 


vV 
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TEST  4  -  MIL.  LOW  CUT  (ZERO  PHASE) 

— filters  •  nm  phase — 

lOM  cut  at  S.  HZi96  n  SLOPE 


REFLECTION  COEFFICIENT 

AT  0B5  IN  LAYER 

1  AT  .0  n  DEPTH 

FREQINH)=  512 

TIH5ECS)=  1000 

niN  =  l.OCHZ 

MAX  =512.00  HZ 

PRIHARIES  ONLY 

NO  AGC.  A5C  = 

1. 

Figure  6.  (Test  4)  Reflections  from  •  one-layer  model 
(flg.  3A)  used  to  display  filter  options.  A,  band  pass 
of  10-50  hz;  B,  band  pass  of  10-500  hz;  C,  D,  E,  F, 
high  and  low  cut  with  minimum  and  zero  phase; 
G,  band  pass  with  6-db  slope  versus  96-db  slope 
used  in  A. 


Test  4.  FUters 


A  single  layer  260  m  thick  is  overlain  by  a 
half  space  wite  the  same  velocity  (1.5  km/sec) 
and  underlain  by  a  half  space  wi^  a  velocity 
of  3.0  km/sec  (fig.  3A).  Ihis  one-layer  model 
was  used  to  test  all  fUter  options.  Layer  and 
half  spaces  have  identical  Q  (200)  and  density 
(1.1).  The  single  reflector  permits  a  study  of 
the  waveform  as  a  variety  of  filter  parameters 
are  changed  (figs.  6A-6G). 

Band  Pass.  Comparison  of  a  minimum- 
phase  band  pass  of  10-50  hz  (fig.  6A)  versus 
10-500  hz  (fig.  6B)  shows  high  frequencies  in 
the  wider  band  pass  (in  the  form  of  a  spike). 
The  dominant  period  on  the  10-50  hz  record 
is  30  msec  (or  33  hz),  about  the  center  of  the 
band  pass.  Band-pass  slope  was  unchanged  in 
figures  6A  and  6B. 

Table  5  is  a  copy  of  the  input  records  used 
to  generate  a  10-50  hz  band-pass  filter  (fig. 
6A).  To  generate  a  10-500  hz  band-pass  ^ter 
(fig.  6B),  one  needs  only  to  change  fl\e  cutoff 
frequency  from  50  to  500. 

Notch  Filter.  If  a  notch  filter  is  desired,  the 
input  is  similar  to  the  band-pass  option  except 
that  the  “filter-type"  is  designated  as  3  for 
both  input  lines.  {See  table  5.) 

High  and  Low  Cut  (zero  phase).  A  high-cut 
filter  (low  pass)  set  at  75  hz  with  96-db  slope 
and  zero  phase  shows  the  reflection  with  a 
dominant  period  of  15  msec  or  66  hz  (fig. 
6C). 

A  low-cut  filter  set  at  5  hz  with  96-db  slope 
and  zero  phase  shows  the  reflection  character¬ 
ized  by  a  sharp  spike  (fig.  6D).  Because  the 
low-cut  filter  essentially  passes  all  .frequencies, 
we  expect  to  see  a  spike.  A  zero-phase  filter  is  ~ 
noncausal;  arrivals  occur  before  the  expected 
arrival  times  of  333  msec  (two-way  distance 
of  500  m  at  1500  m/sec).  Note  that  the 
dominant  arrival  (maximum  amplitude)  oc¬ 
curs  at  333  msec. 

Table  6  is  a  copy  of  the  input  records  used 
to  generate  a  high-cut  filter  (fig.  6C).  To 
generate  a  low-cut  filter  (fig.  6D),  change 
“filter-type”  from  1  to  2  and  change  cutoff 
frequency  firom  75  to  5. 
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Table  5.  Input  data  used  to  test  a  10-50  hz  band-pass  filter  (fig.  6A) 


TEST  4 

-  DIL. 

BAND-PASS 

(10-50  HZ 

EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6. 

CL 

QL 

RHO 

THICKNESS 

1.5 

200. 

1.1 

0. 

1.5 

200. 

1.1 

250. 

3.0 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP -SCALE  MLTPL  AGC  WINDOW 
5  512  1000.  .006  1.0  00  0 
DETECTOR  TYPE  DETECTOR  LOCATKW 
REFLECTION  LAYER  1 

FILTER  PHASE 
1  (MINIMUM  PHASE) 

FILTER-TYPE  E»  SLOPE  CUT-OFF  FREQUENCY 

1  96  50 

2  96  10 


Table  6.  Input  data  used  to  test  a  high-cult  filter  with  zero  phase  (fig.  6C) 

TEST 

4  -  DIL 

.  HIGH 

CUT  (ZERO  PHASE). 

EPS 

SIGMA 

YSIDE  XSIDE 

0.001 

0.1 

40. 

6. 

CL 

QL 

RHO 

THICKNESS 

1.5 

200. 

1.1 

0. 

1.5 

200. 

1.1 

250. 

3.0 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC 

#FREQ  T(MSEC) 

IN/MSEC  AMP -SCALE  MLTPL  AGC  WINDOW 

5 

512 

1000. 

.006  1.0  00  0 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
0  (ZERO  PHASE) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 
1  96  75 


High  and  Low  Cut  (minimum  phase).  The 
model  used  to  study  the  zero-phase  filter 
(above)  was  rerun  with  minimum-phase 
option.  A  high-cut  filter  set  at  75  hz  with 
96-db  slope  and  minimum  phase  (fig.  6E) 
displays  a  slightly  different  waveform  than 
the  zero-phase  filter.  Now  instead  of  the 
amplitude  maximum  at  333  msec,  the  first 
motion  occurs  at  that  time. 

A  low-cut  filter  set  at  5  hz  with  minimum 
phase  (fig.  6F)  also  shows  a  sharp  peak  at  333 


msec.  This  is  consistent  with  the  minimum- 
phase  characteristics  of  high  frequencies 
occurring  at  the  front  of  the  waveform.  Note 
that  the  minimum-phase  filter  passes  a  tail  of 
lower  frequencies  (fig.  6F),  although  the  same 
zero-phase  filter  does  not. 

Table  6,  used  to  generate  the  high-cut  filter 
with  zero  phase,  may  also  be  used  to  generate 
the  minimum-phase  filter.  The  user  only 
needs  to  change  “filter-phase”  from  1  to  0  to 
obtain  figure  6E.  To  obtain  the  low-cut 


I 

p 
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minimum>phaae  seismogram  of  figure  6F,  the 
user  needs  to  also  change  the  filter  type  from 
1  to  2  and  the  cutoff  frequency  from  75  to  6. 

Filter  Attenuation.  Tlie  user  may  select  the 
filter-attenuation  factor  (slope  in  decibel  per 
octave)  for  high-  and  low-cut  filters.  Note  the 
slope  in  decibels  in  a  Butterworth  filter  is 
related  to  the  number  of  poles  (■  DB/6  1). 
Variations  in  slope  were  tested  on  a 
minimum-phase  10-50  hz  band-pass  filter  (fig. 
6A).  With  a  slope  of  96  DB/octave,  the 
10-60  hz  band  pass  is  a  boxcar-type  filter 
yielding  a  reflection  signal  ringing  at  33  hz, 
about  the  centa  frequency  of  the  band  pass. 
But  the  same  filter  with  a  slope  of  6 
DB/octave  yields  a  reflection  with  a  narrow 
spike  (fig.  6G).  The  low  slope  passes  a  large 
range  of  frequencies  to  generate  the  q>ike. 
Figure  6G  is  generated  with  the  input  records 
shown  in  table  5  except  for  a  change  of  slope 
from  96  to  6  db. 

Test  5.  Amplitude  Scale 

The  user  has  available  a  simple  amplitude 
multiplier  ( ASC).  A  two-layer  model  (^.  3B), 
75  m  and  100  m  thick,  was  used  to  test  this 
option.  For  this  model,  reflections  and 
multiples  repeat  at  a  67-msec  interval  For  a 
full-response  option,  an  ASC  ■  100  displays 
six  events  (fig.  5A);  for  an  ASC  *  1,  only 
three  events  are  readily  visible  (fig.  7).  Note 
that  amplitudes  greater  than  0.5  inch  are 
clipped.  Table  4  with  amplitude  scale  changed 
from  100  to  1  will  generate  figure  7. 


TEST  5  -  D2L.  FULL  RESPONSE  (ASC  =  1). 

REFLECTION  COEFFICIENT 

AT  0B5  IN  LAYER  1  AT  .0  fl  DEPTH 

FREQ INN)*  512  T( MSECS )•  1000 

niN  =  l.OtWZ  MAX  *512.00  HZ 

FULL  RESPONSE  (ALL  MULTIPLES! 

NO  AGC.  ASC  *  1. 

NO  FILTERS 


Figure  7.  (Test  5)  Reflectiou  from  a  two-layer  model 
(fig.  3B)  with  smplitude  scale  (ASC)  *  1.  Compare 
with  amplitude  scale  of  100  in  figure  5A. 
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TEST  6  -  D2L.  TEST  OF  AGO. 
RE'^LECTION  COEFFICIENT 
AT  CBS  IN  LAYER  1  AT  .0  i1  DEPTH 
'^REOCNWI=  512  TtMSECSJ=  1000 

HIN  =  l.OOHZ  MAX  =512.00  HZ 
!^0LL  RE5P0NSE  (ALL  MULTIPLES) 

NO  AGC.  ASC  =  1. 

— FIL'ERS  =  MINIMUM  PHASE — 
iIOH  CUT  AT  100.  HZ: 96  OB  SLOPE 


TEST  6  -  D2L.  TEST  OF  AGC. 

reflection  coefficient 

AT  OBS  IN  LAYER  1  AT  .0  M  DEPTH 
FRE(3(NH)=  512  T(MSECS)=  lOOC 

MIN  =  l.OtHZ  MAX  =512.00  HZ 
FULL  RESPONSE  (ALL  MULTIPLES) 

AGC  WINDOH  =  10.  MSECS 

—FILTERS  =  MINIMUM  PHASE— - 
HIGH  CUT  AT  100.  HZ: 96  OB  SLOPE 


a 


;o  1 


:3J  -i 
233  4 
233  M 

!C0  .1 


330  .= 
*O0  .| 
•50  -f 
300  4 
330  4 
300  4 
330  4 

•00  4 

230  4 
300  4 
330  4 
300  4 
330  4 


R1 

R2 

M1 


A 


Figure  8.  (Test  6)  Full  response  from  a  two-layer  model  (flg.  3B)  for  A,  no  AGC;  B,  AGC  with  window  -  10  msec; 
C,  AGC  with  window  ~  60  msec;  D,  AGC  with  window  -  500  msec.  Reflections  R1  and  R2  and  multiple  Ml  are 
identiHed  in  table  4. 
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TEST  6  -  D2L.  TEST  OF  AGO. 
reflection  coefficient 

AT  085  IN  LAYER  I  AT  ,0  H  DEPTH 
FREQ[NW)=  512  Ttf15EC5)«  1000 

MIN  =  l.OIHZ  HAX  =512.00  HZ 
FULL  RESPONSE  (ALL  HULTIPLES) 

AGC  WINOOH  =  60.  MSECS 

—filters  =  MINIMUM  PHASE--- 
HIGH  CUT  AT  100.  M2:96  08  SLOPE 


TEST  6  -  D2L.  TEST  OF  AGC. 
REFLECTION  COEFFICIENT 
AT  085  IN  LAVER  1  AT  .0  M  DEPTH 
FRE01NW)=  512  TIM5ECS)=  1000 

MIN  =  I.OIHZ  MAX  =512.00  HZ 
FULL  RESPONSE  (ALL  MULTIPLES) 

AGC  WINOOH  =  500. MSECS 

—FILTERS  =  MINIMUM  PHASE — 

HIGH  CUT  AT  100.  HZ!96  08  SLOPE 


Flguit  8  —  Continued 
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Test  6.  Automatic  Gain  Control 

Theory.  Normally  we  perform  automatic  suitable  window  of  fixed  length.  Therefore  if 

gain  control  (AGC)  by  dividing  our  time  series  f(t)  is  the  input  signal,  fc(t)  is  the  AGC  signal, 

by  its  own  RMS  value  computed  by  using  a  and  fr(t)  is  the  RMS  signal,  we  have 

fp(t)  =  (/f2(T)  W(t-T)dT)l/2 

and 

to(t)  =  f(t)/t,(t) 


where 

W(t)  =  (  l+OOS(2  TT|t|  /TW)  /  2TW 


(12) 


is  a  window  of  length  TW  designated  by  the  transform  of  f^(t)  and  F(W)  the  Fourier 

user.  For  convenience.  Program  VISP  carries  transform  of  W(t).  F**  is  the  inverse  Fourier 

out  part  of  this  computation  in  the  frequency  transform.  Then 

domain.  Let  F(f^)  denote  the  Fourier 


fc(t)  =  f(t)  /  |F"l[F(f2)F(W)][  1/2 


(  13) 


is  the  algorithm  used. 

Window  Size.  A  two-layer  model  (fig.  3B), 
75  and  100  m  thick,  was  used  to  test  the 
action  of  the  AGC  window  size.  Without  AGC 
and  ASC  *  1,  two  reflectors  (R1  and  R2)  and 
one  multiple  (Ml)  were  observable  (fig.  8A). 
The  time  interval  between  events  is  67  msec. 
A  small  window  (10  msec)  yields  a  noisy 
(ringing)  record  that  does  not  allow  easy 
identification  of  the  early  events  (fig.  8B).  A 
window  comparable  to  the  expected  time 
intervals  (60  msec)  still  has  some  (ringing) 


noise  (fig.  8C).  If  the  window  (520  msec)  is 
one-half  or  more  of  the  time  span  of  the 
entire  record  (1000  msec),  a  seismogram  is 
obtained  similar  to  the  original  record 
without  AGC  (fig.  8D).  The  ENCODE/DE¬ 
CODE  statements  in  subroutine  AGCFREQ 
were  included  to  circumvent  difficulties  with 
significant  figures.  Users  may  have  to  adjust 
this  part  of  the  code  for  their  computers. 

Table  7  is  a  copy  of  the  input  records  used 
to  create  figure  8C. 
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Table  7.  Input  data  used  to  test  automatic  gain  control  (fig.  8C) 


TEST  6 

-  D2L. 

TEST  OF 

AGC. 

EPS 

SIGMA 

YSIDE 

XSII£ 

0.001 

0.1 

40. 

6.0 

CL 

QL 

RHO 

THICKNESS 

1.50 

200. 

1.1 

0. 

1.50 

200. 

1.1 

75. 

3.00 

200. 

1.1 

100. 

1.50 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  *FREQ  T(MSEC) 

IN/MSEC  AMP -SCALE  MLTPL  AGC  WINDOW 

5 

512 

1000. 

.006  1.0  -1  1  60. 

DETECTOR  TYPE 

DETECTOR  LOCATION 

REFLECTION 

LAYER  1 

FILTER 

PHASE 

1  (MINIMUM  PHASE) 

FILTER 

-TYPE 

DB  SLOPE 

CUT-OFF  FREQUENCY 

1 

96 

100 

Test  7.  Computational  Frequency  (NW) 


A  model  with  a  single  layer  (fig.  3A)  was 
used  to  test  the  importance  of  the  user’s 
choice  of  computation  frequencies  (NW)  for 
a  fixed  record  length  of  T  *  1000  msec. 
Figures  9A,  9B,  and  9C  show  signal  quality 
for  computational  frequencies  of  16,  64,  and 
256  hz.  A 10-50  hz  minimum-phase  band-pass 
filter  was  applied  to  the  output.  Comparison 
of  NW  <■  16,  64,  and  256  yielded  seismograms 
that  were  subjectively  evaluated,  respectively. 


as:  not  acceptable,  poor,  and  excellent. 
Although  not  shown  here,  record  quality  for 
NW  *  128  was  good.  For  NW  ■  512  the 
record  was  identical  to  256.  On  a  CDC  Cyber 
170/855  the  times,  excluding  compilation  and 
plotting,  varied  horn  0.088  sec  for  16  hz, 
0.196  sec  for  64  hz,  and  0.621  sec  for  256  hz. 

Table  8  is  a  copy  of  the  input  records  used 
to  create  figure  9. 


Table  8.  Input  data  used  to  test  computational  frequencies  (fig.  9C) 


TEST  7 

-  DIL. 

ROLE  OF 

FREQUENCIES  (NW). 

EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

QL 

RHO 

THICKNESS 

1.5 

2?0. 

1.1 

0. 

1.5 

200. 

1.1 

250. 

3.0 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

«TRC  «FREQ  T(MSEC)  IN/MSEC  AA9-SCALE  MLTPL  AGC  WINDOW 
5  2S6  1000.  .006  1.0  0  0  0 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
1  (MINIMUM  PHASE) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 
1  48  50 
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:EST'  7  -  OIL. 

SOLE  OF  fsequenc:es 

(NWl  test  7  -  OIL. 

ROLE  OF  frequencies  INWI 

TEST  7  -  OIL. 

ROLE  OF  FREQUENCIES  (NM) 

9EfL£CTI0N  coefficient 

pcflectioh  coefficient 

reflection  coefficient 

AT  305  IN  LArCR 

1  AT  .0  H  OEFTH 

AT  005  IN  LATER 

1  AT  .0  N  DEPTH 

AT  005  IN  LAYER 

1  AT  .0  »  depth 

16 

T[HS€C3I>  1000 

FRCQINHI*  64 

nnsECSi- 1000 

FRCQtNNU  256 

T(H5ECSI*  1000 

e  i.o(h/ 

NAX  *  16.00  MZ 

N)N  >  1.0(HZ 

HAX  *64.00  HZ 

niN  --  I.OIHZ 

HAX  *256.00  HZ 

PPIHARICS  ONLY 

PRIMARIES  ONLT 

NC  *GC.  *5C  • 

1. 

NO  AGC.  ASC  * 

1. 

NO  AX.  AX  • 

1. 

phase— 

filters  •  niNimr  ^ase— 

—FILTERS  ■  HINIfW!  PMAX - 

AT  50. 
:ut  at  10. 

HZ; 48  OB  SLOPE 
hZ:40  cm  slope 

HIGH  CUT  AT  so. 
LOW  CUT  AT  10. 

HZi40  00  SLOPE 
hZ:46  00  slope 

HIGH  CUT  AT  SO. 
LOM  CUT  AT  10. 

hZ:48  00  XOPE 

HZ:4B  00  XOPE 

Figure  9.  (Test  7)  Variations  of  computational  frequency  (NW)  for  a  single  reflector  mod'll  (fig.  3 A)  for  A,  NW  ■ 

B,NW-64;C,NW-256, 


Test  8.  Record  Length  (TSEC) 


The  interrelationship  between  the  time 
span  of  the  seismogram  (TSEC)  and  computa¬ 
tional  frequency  (NW)  is  examined  with  a 
six-layer  model  (fig.  3C).  The  record  length 
(TSEC)  was  varied  from  500  to  4000  msec  in 
the  presence  of  a  fixed  computational 
frequency  NW  -  512  hz.  The  minimum 
frequency  is  given  by  1/TSEC,  and  the 
nuutimum  (Nyquist)  frequency  is  given  by 
NW/TSEC.  The  results  ^own  in  figure  10 
demonstrate  that  the  user’s  selection  of  time 
(TSEC)  must  influence  the  appearance  of  the 
seismogram.  For  example,  reflector  R3  of  the 
six-layer  model  (fig.  3C)  is  at  a  depth  of  300 
m  and  has  a  predicted  arrival  time  of  550 
msec.  (See  table  10.)  The  choice  of  record 
length  TSEC  *  500  msec  is  too  short,  and  R3 
aiT'ves  as  a  wraparound  event  at  50  msec  (fig. 


lOA).  Therefore  all  the  arrivals  on  figure  lOA 
other  than  RO,  Rl,  and  R2  are  wraparound 
events. 

A  record  length  of  2000  msec  (fig.  lOB)  is 
sufficiently  long  that  all  events  (RO  to  R6)  are 
observable  at  their  proper  arrival  times. 
Excessive  time  length,  however,  reduces  the 
frequency  of  content  of  the  signal  and  may 
not  be  desirable.  Note  that  a  record  length  of 
TSEC  ■  2000  msec  (fig.  lOB)  gives  an  fmax  of 
only  256  hz  versus  1024  hz  for  the 
seismogram  with  TSEC  ■  500  msec.  One 
should  first  choose  TSEC  to  avoid  wrap¬ 
around  and  then  choose  frequency  (NW)  to 
obtain  the  desired  (maximum)  frequency 
content. 

Table  9  is  a  copy  of  the  input  records  used 
to  generate  figure  lOB. 


TIME  IN  MSEC 
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00 


TEST  8  -  D6L.  TEST  OF  TSEC  (=  5001 
REFLECTION  COEFFICIENT 
AT  OBS  IN  LAYER  1  AT  .0  M  DEPTH 
FREOtNW)=  512  T(nSECS)=500 

HIN  =  2.0CHZ  MAX  =  1024.0(NZ 
PRIMARIES  ONLY 
NO  AGC.  ASC  =  1. 

NO  filters 


TEST  8  -  D6L.  TEST  OF  TSEC  1=  2000) 

REFLECTION  COEFFICIENT 

AT.  OBS  IN  LAYER  1  AT  .0  H  DEPTH 

FREQ(NW)=  512  T(H5ECS)= 2000 

HIN  =  .50  HZ  HAX  =256.00  HZ 

PRIMARIES  ONLY 

NO  AGC.  ASC  =  1. 

NO  FILTERS 


Figure  10.  (Tnt  8)  Variation  of  record  length  (TSEC)  for  a  lix-layer  model  (fig.  3C).  A.  TSEC  -  500  and  B,  TSEC  - 

2000.  Refleetiona  RO  through  R6  are  identified  in  table  10. 
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Table  9.  Input  data  used  to  test  record  length  (fig.  10 A) 


TEST 

8  -  DSL. 

TEST  OF 

TSEC  (=  2000 

EPS 

sicm 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

QL 

RHO 

THICKNESS 

1.5 

10000. 

1.1 

0. 

2.0 

10000. 

1.1 

50. 

3.0 

10000. 

1.1 

250. 

1.5 

10000. 

1.1 

250. 

1.6 

10000. 

1.1 

100. 

3.0 

10000. 

1.1 

250. 

2.8 

10000. 

1.1 

300. 

1.5 

10000. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP -SCALE  MLTPL  AGC  WINDOW 
5  512  2000.  .003  1.0  000 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER -TYPE  I»  SLOPE  CUT-OFF  FREQUENCY 


Test  9.  Arrival  Times  and  Amplitudes 

A  siX'layer  model  with  constant  Q  * 

10,000  and  p  ■  1.1  (fig.  3C)  was  used  to 
examine  arrival  times  and  amplitudes  for  the 
reflection-coefficient  seismogram  (fig.  IIA). 

To  generate  a  noise>free  record,  no  multiples 
were  included.  A  record  length  of  T  «  1200 
msec  avoided  wraparound  problems.  To  avoid 
amplitude  distortion,  the  seismograms  were 
not  filtered. 

Arrival  times  for  the  seven  reflectors  (RO  to 
R6)  were  computed  from  the  velocities  and 


Table  10.  Observed  and  predicted  times  and  amplitudes  of  reflection 
coefficients  for  a  six-layer  model  (fig.  11  A) 


Reflector 

index 

Predicted 

times 

(msec) 

Observed 

times 

(msec) 

Predicted 

amplitudes 

(cm) 

(See  fig.  IIA) 

Observed 

amplitudes 

(cm) 

(normalized) 

RO 

0 

0 

+.14 

.14 

R1 

50 

50 

+.20 

.22 

R2 

217 

215 

-.31 

-.32 

R3 

550 

550 

+.03 

+.03 

R4 

675 

675 

+.25 

+.27 

R5 

842 

840 

-.03 

-.03 

R6 

1056 

1050 

-.22 

-.25 

thicknesses  of  the  layers;  RO  is  from  the  top 
of  layer  1,  and  R6  is  from  the  top  of  the  half 
space.  Theoretical  (predicted)  amplitudes 
were  also  hand  calculated  for  the  model  by 
using  the  usual  relationships  for  reflection 
coefficients  and  two-way  transmission  coeffi¬ 
cients.  Study  of  the  results  (table  10)  shows 
that  the  reflections  arrived  with  correct 
amplitudes,  polarities,  and  times. 

Table  11  is  a  copy  of  the  input  records 
used  to  generate  figure  llA. 


Tint  IN  MSEC 
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TEST  9  -  06L.  ARRIVAL  Tift  TEST. 

RtFLECTlOM  COEFFICIENT 

•T  DBS  IN  LAFEA  I  AT  .0  N  DEPTH 

FOEOiNMi-  siE  TinsEcsi*  mo 

NIN  •  .03  HZ  HAS  •420.57  HZ 

PRIHARIES  tN.T 
NO  AGC.  ASC  ■  I. 

NO  FkTEBS 


TEST  10  -  06L.  arrival  TIME  TEST. 
VERTICAL  OISPLACENENT 
AT  on  IN  LAVER  2  at  .0  M  DEPTH 
FREOIMII*  SI2  TIHSECSI*  1200 

NIN  •  .03  HZ  HAI  •420.67  HZ 

PRiruRIES  ON.* 

NO  AGC.  ASC  •  I. 

NO  filters 


TEST  10  -  06L.  arrival  TIME  TEST. 
PRESSUK 

AT  on  IN  LAVER  2  AT  .0  H  DEPTH 
FREOINNU  SI2  TINXCSI^  1200 

HIN  •  .03  HZ  NAX  •420.67  HZ 

PRIHARIES  OFA-T 
NO  AGC.  ASC  •  I. 

NO  filters 


Figure  11.  (Tests  9  and  10)  Arrival  times  and  amplitudes  for  a  sublayer  model  (fig.  3C)  for  A,  reflection  coefficient; 

B,  vertical  displacement;  C,  pressure. 


Table  11.  Input  data  used  to  test  times,  amplitudes,  and  polarities  of  reflections 

(fig.  IIA) 

TEST  9  -  D6L.  ARRIVAL  TIME  TEST. 


EPS 

SIGM(^ 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

QL 

RHO 

THICKNESS 

1.5 

10000. 

1.1 

0. 

2.0 

10000. 

1.1 

50. 

3.0 

10000. 

1.1 

250. 

1.5 

10000. 

1.1 

250. 

1.6 

10000. 

1.1 

100. 

3.0 

10000. 

1.1 

250. 

2.8 

10000. 

1.1 

300. 

1.5 

10000. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  IFREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  1200.  .003  1.0  0  00 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER-TYPE  DB  SLOPE  COT -OFF  FREQUENCY 
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A; 


Test  10.  Pressure  and  Vertical  Displacement 

The  same  six -layer  model  (fig.  3C)  used  to  top  of  layer  2  (in  effect  at  the  same  depth  as 

test  arrival  times  and  amplitudes  for  the  the  reflection-coefficient  output  but  across 

reflection-coefficient  seismogram  (fig.  IIA)  the  interface). 

was  also  used  to  test  the  waveforms  and  Comparison  of  the  three  outputs  shows 

polarity  of  the  options  to  plot  vertical-dis-  that  vertical-displacement  (fig.  IIB)  and 

placement  and  pressure  seismograms  (figs.  pressure  seismograms  (fig.  IIC)  have  the  same 

IIB  and  llC).  The  algorithm  for  Program  polarity  for  reflectors  R1-R6  but  differ 

VISP  is  so  designed  that  the  vertical-displace-  slightly  in  amplitude.  Reflections  R1-R6  on 

ment  and  pressure  output  may  be  activated  the  reflection-coefficient  seismogram  (fig. 

only  if  the  receiver  is  specified  to  be  within  lOA)  have  the  opposite  polarity.  Because  of 

one  of  the  layers,  not  within  the  half  space.  the  receiver  location,  noted  above,  the  first 

The  reflection -coefficient  plot  is  obtained  arrival  (RO)  on  the  reflection-coefficient 

only  when  the  receiver  is  at  the  bottom  of  the  seismogram  is  a  reflection,  but  the  first  arrival 

upper  half  space.  Therefore,  for  computation-  on  the  other  two  seismograms  is  a  transmitted 

al  purposes,  the  receivers  for  the  three  wave  of  opposite  polarity, 

seismograms  in  figure  11  are  located  as  Table  11  is  a  copy  of  the  input  records 

follows:  the  reflection-coefficient  output  (fig.  used  to  generate  figure  llA.  To  generate 

llA)  is  from  the  receiver  at  the  bottom  of  the  figures  IIB  and  IIC  one  needs  only  to  change 

half  space  (layer  1),  and  vertical  displacement  REFLECTION  to  either  VERTICAL  or 

(fig.  IIB)  and  pressure  (fig.  IIC)  are  from  the  PRESSURE  and  LAYER  1  to  LAYER  2. 

Test  11.  Receiver  Depth 

The  user  may  place  the  receiver  within  any 
of  the  layers.  The  source  is  always  within  the 
half  space,  and  zero  time  begins  when  the 
pulse  encounters  the  top  of  layer  1.  The 
six-layer  model  (fig.  3C),  without  multiples, 
was  used  to  test  this  option.  The  receiver  was 
placed  at  the  top  of  layer  3  (fig.  12A),  and 
vertical  displacement  was  plotted.  The  reflec¬ 
tion-coefficient  option  is  not  available  within 
the  layers  but  only  at  the  bottom  of  the 
upper  half  space.  Table  12  shows  the 
computed  and  observed  times  for  the  primary 
reflections  for  vertical  displacement  (fig. 

12B). 

The  seismogram  (fig.  12B)  shows  the  major 
reflectors  arriving  at  the  computed  times  with 
the  appropriate  polarities.  Table  13  is  a  copy 
of  the  input  records  used  to  generate  figure 
12. 


Table  12.  Arrival  times  for  a  buried  receiver 
(fig.  12) 


Reflections 

Computed 

times 

(msec) 

Observed 

times 

(msec) 

Direct 

108 

no 

R3 

442 

440 

R4 

567 

560 

R5 

733 

730 

R6 

^  947 

950 

s’.V 


1 


«n?t  tmi 
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VI5P;  PLOTS  OF  PARAMETERS 


TEST  11  -  D6L.  RECEIVER  DEPTH. 

CL  lICn/SEC)  QL»  .001  RHO  (GH/CC) 


TEST  11  -  D6L.  RECEIVER  DEPTH. 
VERTICAL  DISPLACEMENT 
AT  DBS  IN  LAYER  4  AT  300.0  M  DEPTH 
FREQ(NW)=  512  T ( MSECS  1=  1200 

MIN  =  .83  HZ  MAX  =426. 67  HZ 

PRIMARIES  only 
NO  AGC.  ASC  =  1. 


Figure  12.  (Te«t  11)  Primary  reflections  for  a  six-layer  model  (flg.  3C)  for  a  receiver  at  the  top  of  layer  4  (where  the 

upper  half  space  is  termed  layer  1). 
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Table  13.  Input  data  used  to  test  receiver  depth  (fig.  12) 


TEST 

11  -  D6L. 

RECEIVER 

DEPTH. 

EPS 

SIGM^ 

YSIOE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

QL 

RHO  THICKNESS 

1.5 

10000. 

1.1 

0. 

2.0 

10000. 

1.1 

50. 

3.0 

10000. 

1.1 

250. 

1.5 

10000. 

1.1 

250. 

1.6 

10000. 

1.1 

100. 

3.0 

10000. 

1.1 

250. 

2.8 

10000. 

1.1 

300. 

1.5 

10000. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  1200.  .005  1.0  000 

DETECTOR  TYPE  DETECTOR  LOCATION 
VERTICAL  LAYER  4 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 


Test  12.  Interpolated  Layers 


It  is  possible  to  have  layer  velocities 
uniformly  increase  (or  decrease)  within 
specified  depth  intervals  without  hand  calcu* 
lating  the  values  (fig.  13).  The  user  enters  a 
sequence  of  values  within  a  sequence  of 


specified  layers:  O,  N,  T  where  O  is  the  flag, 
N  is  the  number  of  layers  to  be  inserted,  and 
T  is  the  total  thickness  of  the  layers.  Table  14 
is  a  copy  of  the  data  used  to  generate  figure 
13. 


Table  14.  Input  data  used  to  interpolate  layers  (fig.  13) 
TEST  12  -  INTERPOLATED  LAYERS. 


EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

RHO 

THICKNESS 

1.5 

200. 

1.1 

0. 

1.5 

200. 

1.1 

75. 

0.0 

5 

500 

1 

3.0 

200. 

1.1 

100. 

1.5 

200. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(IVBEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  1000.  .006  1.0  000 

DETECTOR  TYPE  DETECTOR  LOCATION 
REFLECTION  LAYER  1 

FILTER  PHASE 
1  (MINIMUM  PHASE) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 
1  96  100 
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VISP5  PLOTS  OF  PARAMETERS 

TEST  12  -  OIL.  INTERPOLATED  LAYERS. 


CL  IlCfl/SEC)  QL*  .100  RHO  IGM/CC) 

1  2  3  20  1  2 


1  2  3  20  1  2 


Figure  13.  (Test  12)  Generation  of  interpolated  iayers 
with  uniformly  increasing  velocities. 

Test  13.  Absolute  Depth 

The  user  may  specify  the  absolute  depth 
below  the  bottom  of  the  upper  half  space 
instead  of  specifying  the  intend  thickness  of 
each  layer.  The  absolute  depths  are  entered  as 


negative  values.  Table  15  is  a  copy  of  the 
input  that  in  effect  generates  the  six-layer 
model  (fig.  3C)  used  throughout  this  paper. 
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Table  15.  Input  data  used  to  generate  a  six-layer  model  using  absolute  depths  (fig.  3C) 


TEST 

13-1 

ABSOLUTE 

DEPTH  INPUT 

EPS 

SIGMA 

YSIDE 

XSIDE 

0.001 

0.1 

40. 

6.0 

CL 

RHO 

THICKNESS 

1.5 

10000. 

1.1 

0. 

2.0 

10000. 

1.1 

-50. 

3.0 

10000. 

1.1 

-300. 

1.5 

10000. 

1.1 

-550. 

1.6 

10000. 

1.1 

-650. 

3.0 

10000. 

1.1 

-900. 

2.8 

10000. 

1.1 

-1200 

1.5 

10000. 

1.1 

0. 

9999 

9999 

9999 

9999 

#TRC  #FREQ  T(MSEC)  IN/MSEC  AMP-SCALE  MLTPL  AGC  WINDOW 
5  512  1200.  .005  1.0  000 

DETECTOR  TYPE  DETECTOR  LOCATION 
VERTICAL  LAYER  4 

FILTER  PHASE 
2  (NO  FILTERS) 

FILTER-TYPE  DB  SLOPE  CUT-OFF  FREQUENCY 


INDIANA  GEOLOGICAL  SURVEY  GEOPHYSICAL  COMPUTER  PROGRAMS 

ERRATA 

Geophydeal  Computer  Rrofiim  1  (Occasional  10) 

Pige  9, 19  lines  firom  tbe  bottom  of  the  page: 

Second  line  of  R(MA4)  now  reads 

Second  line  of  R(M.N.4)  should  read  l4-Pa+1^2)4-P(I-i-1^2)+P(M^2)-^P(M^2))/8.0 
hige  9, 4  lines  from  the  bottom  of  the  page: 

Second  line  of  R(M.N.ll)  now  reads  1P(I-20^15)^P(M5^15)4^P(K204+16)-^P(I-*^15^20) 
Second  line  of  R(M^.ll)  ibould  read  lPa-20<)-16)4P(M5^-20)+P(I+20^15)-i-P(K15^20) 

Page  14,  line  6,  which  reads  0(642)^.04007,  may  he  deleted. 

Geophysical  Computer  Program  2  (Occasional  Paper  13) 

Page  11,  line  18: 

Now  reads:  (1,170)1TYPE,Z(I)4U(I) 

Should  read:  (243C)1TYPE,Z(I)4U(0 

Aige  12,  after  line  18: 

hisert:  230  FORMAT  (I1,F4.0,F4.1) 

Geophysical  Computer  Program  3  (Occasional  Paper  14) 

Ibge  12,  line  11: 

Now  reads:  10  A(I+MN)-A(I) 

Should  read:  10  A(M-^K-D-A(N-^K-D 

Geophysical  Computer  Programs  4  and  5  (Occasional  Papers  22  and  23) 

Geophysical  Computer  Programs  4  and  5  require  many  rignificant  figures.  Double  precision  may  be 
needed  on  some  computers.  Indiana  Unirersity  computers  use  60-bit  words. 

Geophysical  Computer  Program  7  (Occasional  Paper  29) 

Subroutine  MYLINE2  has  been  removed  fTom  the  program.  Delete  all  references  to  this  subroutine 
and  read  all  references  to  “11  subroutines"  as  “10  subroutines." 

^e  38: 

Now  reads:  30  X  27  km  region 
Should  read:  31  X  27  km  region 

Ffege39: 

Now  reads:  6  X  40  km  region 
Should  read:  10  X  40  km  region 

Pkge44: 

Now  reads:  distance  200  km 
Should  read:  distance  20  km 

hge  52: 

Now  reads:  as  a  fimction  time 
Should  read:  as  a  fimction  of  time 

Geophysical  Computer  PToparn  9  (Occasional  Fbper  40) 

Ibge  18,  Hne  16: 

Now  reads:  110  THETA-n/2.0 
Should  read:  110  THErAl-n/2.0 


